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ABSTRACT 



Some procedures have been developed to analyze the flowfield of highly 
underexpanded axis ymme trie ring jets operated at high altitudes. The Method 
of Characteristics (MOC) was used to compute the Prandtl-Meyer expansion fan 
and the flow parameters in that region. The MOC may also be used to obtain 
some indications about the repetitive expansion — compression behavior of the 
jet as well as the divergent shape of the expansion part downstream, when 
the ambient pressure goes below certain limits. 

The continuum methods may be used as far as the limit where translational 
equilibrium ceases to exist. The breakdown of the continuum theory may be 
evaluated using the experimental breakdown parameter as proposed by 
G. A. Bird. Beyond this limit, the molecular Direct-Simulation-Monte-Carlo 
method (DSMC) is applied. 

The axisymmetric Monte Carlo Simulation Degins outside the nozzle, where 
the breakdown parameter has a value of approximately 0.05. The actual shape 
of this breakdown limit is a closed lobe surface which starts at the nozzle 
lips, approximately follows a specific Mach wave in the Prandtl-Meyer fan and 
closes back to the axis far downstream. Because our interest is limited to 
the close vicinity of the corners, there the shape of this surface may be 
approximated by a straight line (making a cone in an axisymmetric flow). 

For the simulation purposes , the domain between the breakdown surface and 
the wall is divided into sectors, each sector divided into radial regions and 
each region into simulation cells. Each cell is filled with a number of 
simulated molecules relative to the cell volume and local number density. 
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The simulation is performed for each region separately and contains: 

* molecular motions 

* generation of new molecules to simulate input flows 

* deactivation of molecules to simulate output flows 

* molecular collisions 

Because there is no apriori information about the flow interaction between 
different regions, the whole simulation is done on an iterative basis. A 
first run will supply the output flux from each region. These results are 
used as input data for consecutive runs. 

The program runs as far as the collision frequency is still high and the 
mean free path is low compared with the size of the cells, beyond this limit 
the flow may be regarded as collisionless, and the flux towards the solid wall 
may be computed directly. 
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I. INTRODUCTION 



Gas jets released from spacecrafts and external flows about vehicles at 
high altitudes have a renewed interest in particular with regard to two main 
aspects : 

a. contamination of the spacecraft walls 

b. optical disturbances caused by the plume, 

A spacecraft gas dynamics laser releasing a large quantity of gas is highly 
affected by these two factors and the interest in analyzing them and being 
able to control them have a unique importance in further development, 

A spacecraft laser is assumed to have a long cylindrical shape with the 
optical power output devices installed at one of its bases as shown in Figure 
1 , 

The output gas is released through a ring nozzle, undergoes a fast 
three dimensional axisymmetric expansion, and forms a plume covering the whole 
meridian plane of the vehicle. It widens to large angles of expansion so that 
it may intersect the laser beam. Back scattered molecules may return to the 
wall of the spacecraft causing contamination and degradation of surfaces and 
vehicle parts. 
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Figure 1. The Spacecraft Laser. 

M Q - Mach number at the exit surface ~ 4 the jet gas is composed of two 
species 



heavy molecules 
light molecules 

Altitudes between 200 to 1000 km. 



M G i = 19 

MG2 = 4 
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Continuum flow theory may be used to solve the flowfield and flow 
parameters as far as there is translational equilibrium, it means that 
intermolecular interactions are fast enough to maintain expansion rates. 
Wherever these interactions are too slow, the continuum flow becomes invalid 
and the molecular flow theory should be employed. 

The solution for the continuum regime is computed here by means of the 
Method of Characteristics (HOC) [1,2,3]. The limit where continuum breakdown 
occurs was estimated by the experimental breakdown parameter as proposed by G. 
A. Bird [4]. Beyond this limit, it is proposed to compute the molecular flow 
by means of the Direct Simulation Monte Carlo technique as described in detail 
by Bird [4] . 

For moderate and low pressure ratios (static pressure at the nozzle exit 
to the ambient pressure), an underexpanded jet exhibits a repetitive expansion 
- compression behavior with a geometry depending on the initial Mach angle, on 
the Prandtl-Meyer fan angle, and on the gas specific heats ratio. For lower 
ambient pressure which occurs at higher altitudes, the first compression 
region is pushed out to the envelope of the jet forming the barrel shock. If 
the ambient pressure is low enough, this compression region may disappear due 
to the molecular behavior of the flow [6,11]. 

The breakdown of the continuum theory occurs in a region where the gas 
density and pressure are high compared with the density and pressure in the 
ambient gas. At high altitudes the ratios between these parameters may reach 
10 6 or more. Considering this range of variations, computational 
validity dictates the use of the Direct Simulation Monte Carlo method. In the 
higher density range the jet will be considered as consisting of two species 
of gas, their molecular model will be "the hard sphere molecule" model and 
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ambient gas is not allowed to protrude. In the lower density region the flow 
will be regarded as collisionless. 

In the following chapters we bring the detailed description of the 
computer programs which solve the different parts of the flowfield. 
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II. THE CONTINUUM REGIME 



A. THE TWO DIMENSIONAL ISENTROPIC UNDEREXPANDED JET 

The results brought here are based on the supersonic steady isentropic 
flow theory as described in literature (see for example, Shapiro [1], Liepman 
and Roshko [2], and Owczarek [3]). 

The ranges of paremeters of a jet emerging from a gas dynamics 
spacecraft laser are: 

a. The Mach number at the exit surface Mq » 4. The static pressure 
at the exit plane P Q « 136pa. Ambient pressure (P a mb)> temperature 
(T am b) anc * other thermodynamic properties of amoient gas depend on 
the altitude as shown in Table 1. The jet gas may consist of DF , HF , 
Helium and other species. In the programs we limit the compositon to 
two species: Air and He. (The program allows changes in tne 

composition and types of gas.) 



P 

b. The pressure ratio — — for the minimum required altitude (200 

amb 

km) assures that the jet is highly underexpanded (we show later the 
influence of this ratio on the shape of the jet). 

The following thermodynamic relations are valid as long as the 
compessible flow is isentropic 



= T( 1 


+ ^M 2 ) 


(1) 


= PC 1 


+ JCl M 2 )^ 


(2) 
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(3) 
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where T*p, P-p, and pp* are the total temperature, pressure and density 
(constant for isentropic field). T, P, and p are local temperature, pressure 
and density y is the specific heat ratio of the gas (considered here as 
constant), M is the local Mach number. 

The partial differential equation of motion for supersonic 2-D 
irrotational and isentropic flow is a hyperbolic equation having solutions 
obtained from invariants along characteristic lines. Physical interpretation 
of these lines are the compression or expansion waves which are oriented at 
Mach angles relative to the streamlines. 

Once the directions of the characteristics (waves) are determined 
everywhere in the field, all other parameters may be calculated. 

B. THE TWO DIMENSIONAL PLANAR JET 

Tne compressible supersonic jet flow is characterized by two families 
of characteristics (pressure waves) starting at each corner of the nozzle 
lips. Each of these families of waves forms a Prand tl-heyer fan. The 
streamlines crossing the characteristic waves bend outwards resulting in an 
increase in the flow area. The angle y between the streamline and the 
pressure wave is a function of the local Mach number as 

y = arcs in (1/M) (4) 

The symbol y is the Mach angle. 
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Using the isentropic relations, we can find the relation between the 



turning angle (0) and the local Mach number (M) as 



d9 = - 



M( 1 + M 2 ) 




(5) 
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Integration between the conditions M=1 and M gives the total turning 
angle starting at the throat (where M=l) up to a point with given M. The 
result gives the Prand tl-Meyer function (angle) as 



In the close vicinity of the nozzle lips where the two families of 
characteristics do not intersect with each other there is a "simple region" of 
expansion. There the flow parameters are defined by v and 9 of each 
characteristic line. Further downstream the waves intersect each other. In 
this part of the flow parameters are defined by the two intersecting 
characteristics. A singularity occurs when the initial Mach number is unity 
and a special treatment is required to start the calculations at that point 
(this special treatment has not been brought here). A particular importance 
of the Prand tl-Meyer function is when analyzing the two dimensional flow using 
the hodograph plane. 

C. THE HODOGRAPH PLANE FOR A TWO DIMENSIONAL JET 

The hodograph plane is a representation of the flow parameters in the 
velocity plane. Figure 2 shows a hodograph plane calculated for a simple 
gas (y=1.4). The circles represent constant Mach numbers and constant 




( 6 ) 
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Figure 2. The Hodo graph Plane. 

1-2-3-4-1 defines a jet with Mq> 1 expanding into P a mbl ^ p o> 
resulting simple repetitive expansion-compression. 
l*_ 2 *- 3 *- 3 **- 4 *-i* defines a jet with Mq= 1 expanding into 
P am b2 ^ P o » resulting in a highly underexpanded jet. 
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pressure ratios (“ 5 —) • The epicycloids are the two families of 

x 

characteristics. The angles 9 are the turning angles due to expansion or 

compression (which are both present in supersonic jets). 

To define an isentropic supersonic jet on the hodograph plane it is 

necessary to define the Mach number at the exit plane, the pressure ratios 
p 

(— ) and (— ) , and the specific heat ratio ( y) of the gas, 

T T 

Investigating the shapes of the jet as a function of ranges of parameters 
we may get: 

P amb 

a, simple underexpanded jets for which — is high enough so that 

x 

the two families of characteristics intersect each other. 

b. a critical underexpanded jet for which ' a — is low enough so that 



the intersection between the outer characteristic lines of the two 
familes occur on the outer hodograph circle. 

c. highly underexpanded jets for which a part of characteristics do not 
intersect at all. 

d. expansion into complete vacuum so that there are no reflections from 
the jet boundaries and therefore the compression region disappears. 
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D. THE SHAPES OF TWO DIMENSIONAL JETS 



The following paragraphs further detail the different shapes. 

1 . Simple Underexpanded Jets 

Figure 3 shows the physical plane and the hodograph plane of a simple 

underexpanded jet. If Hq> 1 , an initial turn of the flow 0 O is made within 

the nozzle. An aditional turn of 0 is due to the underexpansion. 0 is found 

by the intersection of characteristics (1-2) with the circle defined by 

^amb/^T • 

When M 0 >1 , the characteristic line 1-2 is described in the physical 
plane by a region in which only one family of expansion waves are present 
(simple region, see transverse line between 1 to 2 in the physical plane). A 
different family of expansion waves forms a second simple region when moving 
between 2 to 3 . At a larger distance from the exit plane, reflected waves 
from the free streamline (jet boundary) cause compression. An ideal 
representation of such a jet is a repetitive pattern of expansion and 
compression. 

For M 0 =l , the tail waves of both families are perpendicular to 
the flow, thus, both lie on line AB . That means that if M 0 =l there is no 
simple region near the exit plane of the jet. The characteristic line 1-2 on 
the hodograph plane becomes a single point A (or B) when located in the 
physical plane. 

2 . Critical Underexpanded Jets 

We define a critical underexpanded jet when point (3) (see Figure 4) 
lies on the limiting circle or P/P^O. This means that there is a core 
within the jet where the pressure approaches zero and Mach approaches 
infinity. This core is theoretically bounded at its upstream side by 
expansion waves and downstream by compression waves. 
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Figure 3 
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Using the theoretical expressions for isentropic ideal flow one 



p p 

may derive the values of — or — as function of M q which causes a jet 



to be critical underexpanded. 



Figure 5 shows results of 



amb 



as functions of M 0 for gases with 



different specific heat ratio. 



3 . Highly Underexpanded Jets 

When Pamb is than the critical values as shown in figure 4, 

point 3 does not exist (there is no intersection between lines 2-3 2 f -3 T ). 

This means that the repetitive reversible expans ion/ compress ion shapes ceases 
to exist. The envelope of the jet starts at an angle defined 'by 9 at the 
nozzle exit plane and approaches an asymptotic angle defined by 9li m (see 
Figure 6). 

In this case we get two (symmetric) groups of characteristics C^-C2 and 
C f |-C f 2 (Figure 6) with no intersection between them. C| and f define the 
inner limit for reflected characteristics, C 2 and C2* define the outer limit. 
As the rest of reflected waves lie between C 2 or C f 2 and the jet boundary, we 
may conclude that a compression region may exist only in a layer along the jet 
boundary. 

because of irreversible effects such as shear stresses, heat transfer due 
to high temperature gradients, or condensation effects in real gases, the 
compression layer may be interpreted as the "barrel shock". 
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Figure 4. The Hodograph Plane for a Critical Underexpanded Jet. 
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Figure 5. Dependence of critical values of ^amb^o 

on the Mach number at the exit plane (M Q ) for different values 
of specific heat ratio (y). 
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c, 




Figure 6. The Hodograph Plane for a Highly Underexpanded Jet 
(Y = 1.4) 

®max “ maximum turning angle 

9 - total turning angle in the specific jet 9 > 9 nia v ./ ? 

9ii m - limiting angle of compression region. 
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Figure (7) shows the hodograph plane for an air jet C Y = 1 • ^ ) expanding 
into an ambient pressure 105 times lower than the static pressure at the exit 
plane. Figure (8) is a schematic description of the jet (the data for 
this jet is given in Table 2). The shape shown in Figure (8) may be compared 
with the barrel shock photograph in page 208 Reference [3]. 

4. EXPANSION INTO A COMPLETE VACUUM - p amb =Q 

This is an extreme situation in which the maximum turn angle of 
streamlines occur near the nozzle exit plane. The theoretical free stream is 
defined only by the theoretical turning angles. All streamlines in the flow 
expanded monotonically towards P»0 without being reflected by the jet 
boundaries . 
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TABLE 2 

CALCULATION OF CHARACTERISTICS ON PHYSICAL PLANE 



State 
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v(M) 
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0+|J 


0-u 


P/P t 
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i 


90° 


90° 


-90° 


0.5283 


2 
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10 


10 


1.435 
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54.15 


-34.15 




3 


0 


20 


20 


20 


1.775 


34.3 


54.3 


-14.3 


.2990 


4 


0 


40 


30 


30 
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27.93 


57.93 


2.07 


.1810 


5 


0 


60 


40 


40 


2.54 


23.18 


63.18 


16.82 


.1035 


6 


0 


80 


50 


50 


3.013 


19.38 


69.38 


30.62 


.55 -1 
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0 


100 


60 


60 


3.595 


16.15 


76.15 


43.85 


.2672 -1 


8 


0 
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67 .7 


67.7 
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53.74 
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20 
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Figure 7. Hodograph piane showing Che points on the mesh of characteristics. 
(Related to physical plan Figure 8). 
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Figure 8. Schematic shape of a highly underexpanded jet. 
(See Figure 7 - The Hodograph Plane) 

The arrows indicate flow direction. 
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E. THE METHOD OF CHARACTERISTICS; COMPUTATION OF PLANAR AND AXISYMMETRIC 
TWO-DIMENSIONAL FLOWS 

For a planar two dimensional flow, the Prandtl-Meyer function v and 
flow direction 9 at any point (3) in the field may be calculated using data of 

two other points (1) and (2) located on characteristic lines that intersect at 

(3) (see Fig. (9)) 

v 3 » J (vj + v 2 ^ + T ^ 9 1 “ 9 2^ ( 7 ) 

9 3 = I (v l " V + T (9 1 + 9 2 } (8) 




Figure 9. The calculation of v and 9 for point (3) is based on data for 
points 1 and 2. 
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For axisymmetric , two-dimensional flow, Liepmann and Roshko [ 2 J developed 
expressions for finding the Prandtl-Meyer function (v) and the flow direction 
( 8) which are given by: 



sin0 



-T ( h 



V 



+ T< 9 i - 



V 



+ ~ [sin 



1 



sin0„ 



+ sin 



An 23 J 



(9) 



i (v i 



sin9 



V 



(9 i + V + i [sin 



1 



sin0 2 

A?i 3 - sin w 2 — - An 23 J 



( 10 ) 



The angles and subscripts are shown in Figure (10), 




Figure 10. Calculation of 0 and v for axisymmetric flow. 
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It is obvious that for the axially symmetric flows the increase in the 
radius causes the increase in the flow cross section and influences the flow 
direction and the Prandtl-Meyer function* These facts have been taken into 
account when developing the "compatibili ty equations” (9,10). 

Using equations (9,10) we have developed a computer program which 
enables the calculation of the jet flow for a two dimensional and for an 
axially symmetric geometry (ring jet). 

The listing of the program, the program description and some results are 
given in Appendix A. 
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III. THE BREAKDOWN OF THE CONTINUUM THEORY 



A. GENERAL CRITERIA 

As described in detail by Bird (Chapter 1 Reference [4]), the validity 
of the continuum approach has been identified with the validity of the Navier 
Stokes equations. This requires that the Knudsen number K n = X/L should be 
small compared with unity (X is the mean free path and L is a scale length for 
the specific flow field). For K n larger than a certain limit (between 0.01 
to 0.1 depending on the required accuracy) a microscopic approach is 
necessary. 

For small values of L , the microscopic approach may lead to statistical 
fluctuations of the results due to the small number of molecules participating 
in the flow processes. In figure (11) which was reproduced from Bird’s book 
[4, figure 1.6], the regimes of rarefied flow and high fluctuations are 
depicted. The flow around the jet-air boundaries near a spacecraft is 
generally rarefied (high Knudsen number), but has insignificant fluctuations. 
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Figure 11. Limits for continuum approach and microscopic approach 
(d-3.7 x 10-10 m ). 
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B. THE EMPIRICAL CRITERION 



The Method of Characteristics (HOC) was used to compute the jet flow 
and the results obtained from the computer program "AXSYM" are valid as long 
as the continuum flow theory is valid. 

The continuum flow requires that the mean free path should be negligibly 
small in comparison with the scale length of the macroscopic flow variations. 
The classical theory for Prand tl-Meyer expansion may therefore be expected to 
fail at progressively larger distances from the nozzle lips as the gas density 
decreases with the increasing flow angle and Mach number. The empirical 
criterion for the breakdown of continuum flow in steady expansion flow [4] is 
that 



P 




d^ 

dS 



0.05 



( 11 ) 



where 



q = stream velocity 
p = density 

v = molecular collision frequency 




absolute change in density while moving a distance 



dS along a streamline 

Introducing the breakdown parameter P into the program, gives the 
definition of the boundary where the flow should be calculated oy means of the 
molecular flow theory, i.e., by solving the Boltzmann equation. 
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For an underexpanded jet with a high initial Mach number, 
the breakdown surface is nearly a streamline. Furthermore, the range of flow 
parameters for the present problem are such that the simple region extends to 
very large distances and near the nozzle lip the breakdown limit may be 
approximated by a straight line. 

For the axisymmetric jet there is no simple region, however, for the 
region of interest it may be regarded as linear. 

The method proposed in the present work for solving the flow behind the 
breakdown boundary is the Direct Simulation Monte Carlo (DSMC). For tnis 
purpose, a computer program "SIMUL" was developed. In the following chapter 
we describe the algorithms required for the specific problem, the geometry and 
the data organization. Detailed program description is given in Appendix (B) . 



i 
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IV. THE MOLECULAR FLOW IN AN AXISYMMETRIC RING JET 



A. GENERAL CONSIDERATIONS 

The part of the field in which the jet may be calculated by means of the 
continuum theory was described in Chapter II. There we calculated also the 
boundaries where continuum theory becomes invalid and molecular calculation 
should be employed. In fact, the molecular theory and the molecular Boltzmann 
equations are universal and hold for the entire flowfield. However, 
computational requirements make the Boltzmann equation impractical for the 
upstream flows. Therefore we limit our solution only to the part of the flow 
beyond the region where continuum breakdown occurs. 

As a result obtained from MOC solution the "breakdown" , i.e., the locus 
where the breakdown parameter p has values between 0.03 to 0.06, for the 
region close to the nozzle lips this boundary may be approximated by a 
straight line (for axis ymme trie flow this line is the envelope of a cone, see 
Figure (12)). 

For the specific jet and gas, the breakdown occurs in a region where the 
number density is in a range of 10^1 molecules/m 3 . For ambient gas at an 
altitude of 200 km the number density is 10 15 and decreases to a range of 10 1 * 
at 1000 km. In order to be able to express this vast change in a simulation, 
we would need to have an extremely large number of cells which would be 
impossible to store in a computer. To overcome this problem we are required 
to make a less exact formulation which enables the production of results, 
having to pay the penalty of "smearing" the steep gradients and obtaining 
averages within layers of simulated cells. Unfortunately it is impossible to 
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predict how far the simulated results will be from the exact solutions. These 
comparisons have to be made after getting final results of this simulation. 

1 . The Direct Simulation Monte Carlo Method 

The direct simulation Monte Carlo Method is a technique for a computer 
modeling of a real gas by some thousands of simulated molecules. The velocity 
components and the coordinates of the simulated molecules are stored in the 
computer and are modified with time as a result of collisions and boundary 
interactions. A detailed description of some problems and their solutions 
by means of direct simulation is given in [4]. 

To follow the molecular motion it is necessary to divide the simulated 
domain into a network of cells . The size of _a cell mus t be such that the 
change in flow properties across each cell is small . The time is advanced in 
discrete steps DTM, such that DTM is small compared with the mean collision 
time per molecule . If there is a flow going through the domain, DTM should De 
small compared with the mean time required for the mean flow to cross trie 
cells.* Both cell size ( ( DR) , ( DDALFA) - radial size and angular size as they 
appear in the program) and DTM may vary in the simulation with position and 
time . 

Applications such as free jet expansion in which large gradients of flow 
properties are expected, may require a very large number of cells for the 
simulation. In these cases the computer memory requirements to store cells T 
data and molecules 1 data may exceed the available computer storage. A 



*If DTM is chosen to be very small compared with the mean time between 
collisions then the simulation will require a very large number of runs such 
that the number of collisions will be sufficient. If DTM is large the 
molecules are washed out by the mean flux and there is no time for the 
collisions to influence the flow. 
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solution for this problem is to divide the simulation space into smaller 
regions and to run the simulation for each region separately. If there is an 
interaction between different regions, which apriori is undefined, the 
solution should be found iteratively. (That means that each run will provide 
data for consecutive runs and the procedure should be repeated until the 
results converge to a steady solution. 

The computation of a representative set of collisions based on mean 
collision time per molecule is invalid for a computerized simulation because 
of the large computer time and computer memory requirements. Instead, the 
method proposed by Derzko which is described in details by Bird [4] may be 
employed. Following this method, an averaged mean time between collisions of 
species L with species M for a cell is calculated. The number of 
collisions of each type (L.M species) is such that the collision time 
counters are kept concurrent with the overall time parameter. The L.M 
collision time for a cell containing and Njj molecules with collision 

cross section a^i > number densities n^ and njq and relative velocity C r , 
is given by 



where LP 
for the L 



_ LP 1 

At - — 

C l L °LM n M C r 



+ 



HP 



N M °LM n L C r 



( 12 ) 



and MP are the probabilities that the collision will be effective 
and 11 molecules respectively. 



31 



B. THE GEOMETRY OF THE SIMULATED DOMAIN, SECTORS, REGIONS AND CELLS 

Figure (13) shows a cross section of the simulated domain for the 
axisymmetric (ring) flow. Points A and A 1 are the nozzle lips. Starting at 
"A" and assuming the "breakdown” boundary to be a straight line, we obtain the 
cross section of the molecular domain as a sector defined by LAM. The solid 
wall is defined by AL. The arc LM may be assumed to be far enough so that the 
pressure along it may be assumed to equal the ambient pressure. Molecules 
originated in the jet cross the breakdown boundary with a velocity, direction, 
temperature (and other thermodynamic properties) as found from the continuum 
solution. 

The molecular domain LAM is divided into secondary sectors, and each of 
these are divided into several radial regions making the "simulation 
regions" . 

Because we have no apriori information on how the expansion occurs, the 
angle of each sector (which mainly is in the direction of the expansion 
gradients) is left to be a result of the internal calculation. 

Each region is divided into NRD radial divisions and NAD angular 
divisions making a network of NAD*NRD simulation cells. The angle DALFA of 
all regions in a sector is constant. Taking NAD constant for all regions in 
a sector, we get the angle of a cell DDALFA constant. Defining the radial 
size of a cell DR as constant we get a cell cross section area proportional 
to the radius R measured from the nozzle lip (point A). 
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Figure 13. Cross section of the simulation domain, definition of sectors, 
regions, cells and coordinates. 
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The size of a cell: In order to get accurate simulated results it is 

recommended to define the size of a cell (DR and R*DDALFA) small compared with 
the mean free path of the molecules X (typical DR=X/3). However, as we do 
not expect to get large changes in flow parameters along the radius we may 
allow DR be much larger than X/3. The angular size of the largest cell in a 
sector should comply with this requirement, but because of the computer 
limitation it is set to be equal to 5*X. This will be the basis for defining 
DALFA for each sector, 

C. INITIAL NUMBER OF MOLECULES IN CELLS 

A '‘reasonable" number of molecules in a simulation is several thousands 
(a larger number, which is better, may be used for simple problems or when 
using a single user computer with large user memory space). The initial 
setting of molecules in cells is usually based on a guess of the number 
density in the specific cell. (The number of molecules in cells will change 
during the simulation according to the input/output calculated fluxes to the 
specific region). 

The number density and the size of a cell are specified only in three 
dimensional flows. When applied to a two dimensional flow the simulation may 
be regarded as applying to an arbitrary thin slice of the real flow. In the 
axis ymme trie flow we define the width of a cell by the angle DFI as shown in 
Figure (14), constant within a region. 

The initial number density in cells of a given region is set constant. 
Defining the total number of simulated molecules in the region the number of 
molecules in each specific cell becomes a function of DFI. For example, 
assume that we limit the number of molecules in the smallest cell in the 
region to 15 then: 



34 




Figure 14. Variation of the angle DFI. 

To maintain the number of simulated molecules within computational 
limits, the 'width' of each region defined by DFI is such that: 

MIN * number of molecules in smallest cell in a region 
MIN =* VOLUME (smallest cell) * number density 

« f (R, ALFA, DALFA) * DFI * number density for flux 
calculations 

DFI is a weighting factor. 
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DFI*(contant) * (number density) = 15 
Other cells contain the initial number of molecules proportional to their 
volumes • 

D. DEFINITION OF INPUT AND OUTPUT FLOV/S FOR A REGION 

The cross sections of all regions (except those regions near the nozzle 
lips) are quadrilateral. Through the sides of the region molecules are 
allowed to enter or to leave according with the boundary conditions or as a 
result of molecular velocity. For the first sector, near the breakdown 
boundary the input flow (FWP1 and FWP2 see Figure (15)) is defined by the 
results from the continuum flow. FEN1 , FNN1 , FSN1 , etc. are results of 
counting and averaging the outgoing molecules (the different vector names will 
be explained in Appendix B) . 

For the neighbor regions these output fluxes become inputs and have to 
be adjusted according to the differences in the angle DF1 of the different 
regions . 

The simulation starts with regions in the sector near the 
breakdown boundary. At this time there is no data for input flows through 
faces E and N of the cell. An additional run of the whole program is required 
in order to take these calculated flows into account. If the accuracy of the 
results is important we may run this type of iteration several times until 
the results become stable. (Only after running the program for the whole 
domain once we shall be able to evaluate the importance of these iterations.) 
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Figure 15. Definition of input and output flows of species 1 to a region (KR) 
in a sector (KS). 

(For species 2 the flux names will change as follows: 
instead FWN1( ) * FWN2( ) 
instead F0W1( ) > F0W2( ) 
etc. 
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E. COLLISIONLESS FLOW 



In several sectors near the breakdown boundary we may find a high number 
density and the mean free path small compared with the size of a cell* There 
the calculated collisions are expected to have an influence on the flow 
parameters. For wider expansion angles the collisions become rare mainly 
because of the decrease in the density. In the ambient gas the mean free path 
(for 200 km altitude) is 240 m. Comparing this number with the size of the 
simulated domain may lead to the conclusion that there the flow may be 
regarded as collisionless. 

We may define a limiting line in the flow where the collisions become 
insignificant. Consequently, molecules crossing this limit will in fact 
continue moving in straight lines; a part of them reach the solid wall. 

Introducing this idea of the collisionless flow we may reduce the 
computation time and the memory requirements. 

F. TWO DIMENSIONAL PLANAR FLOW VS. AXISYMMETRIC FLOW 

The cell dimensions are completely specified only in three dimensional 
flow. When applied to two dimensional flow, the simulation may be regarded as 
applying to an arbitrarily thin slice of the real flow. The thickness of the 
slice may be chosen such that the number of simulated molecules complies with 
the cell volume and the physical number density. For the axisymmetric flow we 
have defined the angle DFI as the third coordinate so that the volume of the 
cell is completely specified. 

Once the geometry is defined, the simulation may be accomplished and 
there is no difference if doing it for two dimensional or for axisymmetric 
f lows . 
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APPENDIX A - THE AXSYM PROGRAM 



A. 1 DIFFERENT REGIONS IN THE JET 

For the two dimensional jet with initial Mach number greater than unity 
the different regions are shown in Figure (16). 



Figure 16. The three regions in an underexpanded jet. 

For planar 2-D flow region 1 is a uniform flow core, region 2 is a simple 
region in which only one family of characteristics define the flow, and region 
3 which contains the intersection of the two families of characteristics. 
Because our intention is to find solutions for highly underexpanded jets with 
very low ambient pressure, the calculation of further downstream flow is not 
necessary. 

For the axisymmetric ring jet, we use the same definition for the 
different regions however, in this case none of the three regions has uniform 
flow and is not a simple region. 
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As shown in Equations (9,10) the PM function ( v) and the flow direction 
( 9) of a point at location I,j may be calculated from the v and 9 of two 
upstream points (I-l,j) and (I,j-1). Later, from the PM function at the new 
point we may derive the local Mach number, the local pressure, temperature, 
velocity and other thermodynamic parameters as required. 

Definition of the mesh of points for the different regions is shown in 
Figure (17). 

A . 2 PROG RAM FLOWCHART 

A simplified flowchart for the MOC program is shown in Figure (18). The 
program is designed to solve both axisymmetric as well as two dimensional flow 
for kD = 2 it solves two dimensional flow 

for kD = 3 it solves axisymmetric flow (This is also the default 
condition) 

Initial data such as Mach number and pressure at the exit surface, ambient 
pressure and jet gas parameters are input data. 

Output data contains the following for each mesh point: 

Mach number, coordinates of mesh point (R,X), flow direction (TETA) , 
pressure, temperature, local velocity, Knudsen number based on the 
distance between two points along a streamline, mean free path and 
breakdown parameter as defined by Bird.* 

For each of the three regions, we start with precalculated boundary 
conditions enabling the calculation of the ilach angles, coordinates of mesh 
points and distances d£ and dq as described in [2]. 



*For the exit plane instead the Bird's breakdown parameter, we calculate the 
ratio between time per three collisions and time of motion. Sometimes this 
ratio may be regarded as a measure of the breakdown of the continuum theory. 
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Figure 17. Indexing of mesh points for the different regions in the 'AXSYM' 
program. 
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Figure 18. AXSYM Program Flowchart 
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The number of characteristics used in the program is arbitrary and 
depends on the required resolution (it may affect also the accuracy of the 
results and the amount of computation). We start with 20 characteristics 
along the exit cross section and with 50 characteristics in the Prand tl-Meyer 
fan thus a total of 70 characteristics of each family are calculated. In 
region 1 there are 20 left running and 20 right running characteristics. In 

region 2 there are 20 right running and 50 left running characteristics. In 

region 3 there are 50 left running and 50 right running characteristics. 

In region 3 we limit the calculation where the two characteristics 

defining a new mesh point intersect at an angle smaller than the 

computational accuracy. In fact this occurs far downstream where continuum 
theory becomes invalid. 

After defining the mesh geometry (successively) we calculate the 
Pranatl-Meyer function and flow direction, using equations (9,10). 

The local Mach number is an implicit function of the Pradtl-Meyer angle. 
It is calculated by iterations with an initial guess of local Mach number 
set equal to a precalculated Mach number at an adjacent point, and the slope 
of the function as given by Equation (5). Figure (19) shows the iterative 
procedure for evaluating the Mach number at each mesh point. 
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Figure 19. Iterative procedure for liach number calculation. 

Once the local Mach number has been found all other flow 
parameters may be defined using the Equation (1,2,3). 
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Figure 20. The mesh of characteristics in Region 2 and Region 3 
Axi symmetric ring jet. M 0 = 4. Altitude=200km. 
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Figure 21. The mesh characteristics in Region 2 and Region 3 
Axisymmetric ring jet. M Q = 2. Altitude = 200km. 
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A. 3 PROGRAM 'AXSYM' LISTING 



AXS00010 
AXS00020 
AXS00030 

THIS PROGRAM CALCULATES THE ISENTROPIC EXPANSION OF A JET BY MEANS OFAXS00040 
THE METHOD OF CHARACTERISTICS. 

FOR A TWO DIMENSIONAL JET 'KD' SHOULD BE SET EQUAL TO 2 

FOR AN AXISYMMETRIC RING JET 'KD' SHOULD BE SET EQUAL TO 3 



$ JOB 

C PROGRAM AXSYM 
C 
C 
C 
jC 
c 
iC 



IMPLICIT REALX8(A-H,0-Z,$) 

DIMENSION TETA(20,50), AM( 20,50 ),R( 20,50 ),X( 20,50 ),PM( 20, 50) 



AXS00050 
AXS00060 
AXS0007 0 
AXS00080 
AXS00090 
AXS00100 



DIMENSION AMCOR ( 20 ,20),TETAC(20,20),XC(20,20),RC(20,20), PMC( 20 , 20)AXS00110 



DIMENSION DENSF( 20 , 50 ) 

DIMENSION AMX( 50 , 50),TETAX(50,50),XX(50,50),RX(50,50), PMX( 50 , 50) 



AXS00120 
AXS00130 
AXS00140 
AXS00150 
AXS0016 0 
AXS0017 0 
AXS00180 
AXS00190 



TETA IS THE FLOW ANGLE (RADIANS) MEASURED FROM X AXIS. 

AM IS THE MACH NUMBER 

R IS THE RADIUS (NORMAL TO THE WALL) 

X IS THE AXIAL LOCATION (PARALLEL TO THE WALL) 

PM IS THE PRANDTL MEYER FUNCTION 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS00200 
C THE FOLLOWING IS DATA FOR THE SPECIFIC PROBLEM AXS00210 

PAMB = 8.4736E-5 AXS00220 

KD = 3 AXS00230 

C FOR TWO DIMENSIONAL FLOW KD = 2 ,-FOR AXISYMMETRICAL FLOW KD=3 AXS00240 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS00250 



c 










AXS00260 


c 


CONSTANTS 








AXS00270 




PI = 3.141593 








AXS00280 




BOLTZ = 1.38032E-23 








AXS00290 




AVOG = 6 . 0225E+26 








AXS00300 




RG = 8314.3 








AXS00310 


c 


EXIT SURFACE 








AXS00320 




AMO = 4.00 








AXS00330 




TO = 300.0 








AXS00340 




PO = 136.0 








AXS00350 


c 










AXS00360 




R1 = 2.5 








AXS0037 0 




XI =0.5 








AXS00380 


c 










AXS0039 0 


c 


GAS DATA 








AXS00400 




GAMA = 1.535 








AXS00A10 




DIAM = 2.95E-10 








AXS0042Q 




GM =17.0 








AXS00430 




RJ = RG/GM 








AXS00AA0 




CXS = PIXDIAMXDIAM 








AXS00A50 




GMM = GM/AVOG 








AXS0046 0 


c 










AXS00470 


c 


MESH DEFINITION 








AX500A80 


c 


N1 CHARACTERISTICS 


FROM 


THE 


EXIT PLANE 


AXS00490 


c 


N2=CHARACT ERISTICS 


FROM 


THE 


CORNERS 


AXS00500 




N1 = 20 








AXS00510 




N2 = 50 








AXS00520 


c 










AXS00530 


c 


CONSTANTS FOR THE ISENTROPIC EXPANSION 






AXS00540 




A1 = (GAMA-1 .0)/GAMA 








AXS00550 




B1 = 1.0/A1 








AXS00560 


c 










AXS00570 




A2 =DSQRT( (GAMA+1 . )/(GAMA-l . )) 








AXS00580 




B2 = 1.0/A2 








AXS00590 


c 










AXS00600 




A3 = (GAMA-1. )/2. 








AXS00610 


c 










AXS00620 


c 


C = MACHXMACHXA3 + 1. 








AXS00630 


c 


D = MACHXMACH -1. 








AXS00640 



CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS00650 
Cxxxxxxxxxxxxxxxxxxxxxxxxxxxx*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS00660 
C AXS00670 

AXS00680 
' AXS00690 



C DEFINE STAGNATION PARAMETERS 

c 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx AXS007 00 



C = (1.0+A3XAM0XAM0) 
C PRESSURE 



AXS00710 

AXS00720 



PN = CxxBl 
PSTG = POXPN 
C TEMPERATURE 
i TSTG = TOXC 

C DENSITY 

: DO = PO/CRJXTO) 

DSTG = PSTG/CRJxTSTG) 

C 

Cxxxxxxxxxxxxx 

WRITEC6,1)PSTG,TSTG,DSTG,DIAM,GM 

1 FORMATC 'O', 'STAGNATION PRESSURE 2 ', E12 . 5, 
1* DENSITY 2 ', E12. 5, ' MOL . DIAM= ' , E12 . 5, ' 

WRITE C6,2)P0,T0,D0,AM0 

2 FORMATC 'O', 'EXIT PLANE PRESSURE 2 ' , E12.5, 
1' DENSITY=',E12.5, ' MACH= ' , FI 0 . 5 ) 



TEMPERATURE 2 ' 
MOL . MASS 2 ' , FIO 

TEMPERATURE 2 ' 



AXS0073C 
AXS0074C 
AXS0075C 
AXS0076G 
AXS0077G 
AXS0078G 
AXS0079C 
AXS0080C 
AXS0081C 
AXS0082C 
FIO . 5, AXS0083C 
5) AXS0084C 
AXS0085C 
FIO . 5, AXS0086 C 
AXS0087 C 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0088C 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0089C 



C AXS8090C 
C DEFINE FREE STREAM PARAMETERS AT THE RIGHT CORNER AXS0091C 
C AXS0092C 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS0093C 



C MACH NUMBER 

FSM 2 DSQRTCCCPSTG/PAMB)xxAl-l.)/A3) 

C TEMPERATURE 

FST 2 TSTG/C1.+A3XFSMXFSM) 

C DENSITY 

FSD 2 PAMB/CRJXFST) 

C MACH ANGLES FOR HEAD AND TAIL OF FAN 
AMIT =DARSIN(1 ./AMO) 

AMIH =DARSIN(1 ./FSM) 

C PRANDTL MEYER FUNCTION FOR HEAD AND TAIL OF FAN 
D1 =DSQRT(AMO*AMO-l . ) 

D2 =DSQRT ( FSMXFSM-1 . ) 

PMH 2 A2*DATAN(B2XD2)-DATAN(D2) 

PMT 2 A2XDATAN(B2*D1)-DATAN(D1) 

C EXTERNAL TURNING ANGLE (FREE STREAM ANGLE)=EXTA 
EXTA 2 PMH - PMT 
C PRANDTL MEYER FAN ANGLE PMFA 
PMFA 2 EXTA - AMIH + AMIT 
C 

C x x x x x xx x x x * x * * 

C WRITEC6 , 3)PAMB , FST, FSD, FSM 

C 3 FORMATC 'O', 'FREE STREAM PRESSURE 2 ', E12 . 5, ' 
C 1' DENSITY 2 ' ,E12. 5, ' MACH 2 • , FI 0 . 5/// ) 



AXS0094C 
AXS0095C 
AXS0096 C 
AXS0097C 
AXS0098C 
AXS0099C 
AXSOIOOC 
AXS0101C 
AXS0102C 
AXS0103C 
AXS010AC 
AXS0105C 
AXS0106C 
AXS0107 C 
AXS0108C 
AXS0109C 
AXS0110C 
AXS0111C 
AXS0112C 
AXS0113C 
AXS0114C 

TEMPERATURE 2 ' , FIO .5, AXS0115C 

AXS0116C 



Cx*xxxxxxx*xx**x*x*x*xx*x*x*x**x*xx*x*x**x****xx*x****x***xxxxxxxxx*xx**AXS0117C 



Cxxx*x*x*xxxx*x*xx*xxxx*x*xxxxxxxx*xx*xx*x*xxxxx*x*xxxx**x**x*x**x*x*xx*AXS 0118 C 



C AXS0119C 
C DEFINE THE MESH OF CHARACTERISTICS-AND FLOW PARAMETERS AXS0120C 
C AXS0121C 



CXXXXXXXXXXXXXX*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0122G 



C THE CORNER POINT - LEFT RUNNING CHAR. 

C PMFA IS DEVIDED INTO N2 C=50) NONEQUAL DIVISIONS 
C 

RAT 2 CFSM/AM0-l)/0.02 
RAT2 2 FLOAT C N2-1 ) 

El 2 DLOGC RAT)/DL0GCRAT2) 

C 

WRITEC 6,4) 

4 FORMATC '1' , 'PRANDTL MEYER FAN LINE MACH 

1 P.M. ANGLE TETA'/) 

DO 20 N=2.N2 
EN2 2 FLOAT C N-2 ) 

EN1 2 FLOAT C N~1 ) 

AMB 2 AM0xCl.+.02xCEN2)xxEl) 

AMF 2 AM0XC1. + .02XCENDXXE1) 

DELM 2 AMF-AMB 
AMC1,N) 2 AMF 

D1 =DSQRTCAMC1,N)XAMC1,N)-1 . ) 

PMC 1 , N) 2 A2XDATANCB2XDD-DATANCD1) 

AMI =DARSINC1./AMC1,N)) 

TETACl.N) 2 PI/2.-CPMC1, N)-PMT) 

C TETA IS THE FLOW ANGLE ON THE CHARACTERISTICS AT THE CORNER 



AXS0123C 
AXS0124C 
AXS0125C 
AXS0126D 
AXS0127C 
AXS0128G 
AXS01290 
AXS01300 
AXSO1310 
AXSOI320 
AXS01330 
AXS01340 
AXS0135C 
AXS0136 0 
AXS01370 
AXS01380 
AXS01390 
AXS01400 
AXS01410 
AXS01420 
AXS01430 
AXS01440 



48 



R(1,N) = R1 
X(1,N) = XI 
ALFAL = TETA( 1 , N)+AMI 
ALFAL IS THE ANGLE OF THE 



LEFT RUNNING CHARACTERISTICS 



PRESSURE, TEMPERATURE AND DENSITY VARIATION AT THE CORNER 
C = AM(1,N)XAM(1,N)XA3+1. 

PRES = PSTG/CCxxBl) 

TEMP = TSTG/C 

DENSF( 1 , N) = PRES/(RJXTEMP) 



Cxxxxxxxx 

WRITE (6,5) 
5 FORMAT ( » 

CXXXXXXXX 

20 CONTINUE 
C 



N,AM(1,N),PM(1,N),TETA(1,N) 



' , 115, 3E20 .5) 



AXS01450 
AXS01A60 
AXS01470 
AXS01480 
AXS01490 
AXS01500 
AXS01510 
AXS01520 
AXS01530 
AXS01540 
AXS01550 
AXS01560 
AXS0157 0 
AXS01580 
AXS01590 
AXS01600 
AXS01610 



CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxaXS01620 

C AXSD1630 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS0l640 
WRITE(6 ,10) AXS01650 

10 FORMAT ('O', ' I J MACH R X TETA AXS01660 

1TEMP PRESS VELOCITY KNUDSEN MFP ND AXS01670 

1 P’) AXS01680 

C AXS01690 

C AXS01700 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS01710 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS01720 

C AXS01730 

C DEFINE THE PARAMETERS AT THE EXIT SURFACE (PLANE) AXS01740 

C AXS01750 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS01760 



c 

c 



25 

THE 



c 

c. 

c 

c 

c 

c 



DO 25 J = 1,N1 
AMCOR( 1 , J ) = AMO 
D = DSQRT(AM0XAM0-1 . ) 

PMC( 1 , J ) = A2XDATAN(B2XD)-DATAN(D) 

TETAC( 1 , J ) = PI/2. 

EXIT PLANE IS DIVIDED INTO (Nl-1) DIVISIONS (N1 POINTS) 
DO 30 J = 1 , 20 

SOUND = DSQRT (GAMAXRJXTO) 

VELO = AMCOR(l, J)XSOUND 
DNO = PO/(BOLTZXTO) 

FPO = . 7 07/ ( DNOxCXS ) 

DX = X1X2 ./FLOAT ( N1 ) 

AKNO = FPO/DX 

XC(1,J) = Xlx( 1 . -FLOAT ( J-l )/ FLOAT (Nl-l)x2. ) 

CENTR = DABS(XC( 1 , J ) ) 

IF(CENTR.LT. 0.001) XC(1,J) = 0. 

RC( 1 , J ) = R1 



PARAMETER IS 
AND THE. FLOW 



KZ=lxxxxxxxxxxxxxxxxxxx 
IF (J.GT.l) GO TO 29 
AT THE EXIT PLANE THE BREAKDOWN 
THE RATIO OF THE COLLISION TIME 
THE MESH DIMENSION. 

SCALE = DXXDTAN(DARSIN(1 ./AMO) )/2. 

TIME1 = SCALE/VELO 

TIME2 = 3./(<4.xCXSxDNOxDSQRT(BOLTZxTO/(PIxGMM))) 
P =TIME2/TIME1 
DENSF (1,1) = DO 
29 CONTINUE 



EVALUATED BY MEANS OF 
TIME. LENGTH SCALE IS 



AXS01770 
AXS01780 
AXS01790 
AXS01800 
AXS01810 
AXS01820 
AXS01830 
AXS01840 
AXS01850 
AXS0186 0 
AXS0187 0 
AXS01880 
AXS01890 
AXS01900 
AXS01910 
AXS01920 
AXS01930 
AXS01940 
AXS01950 
AXS01960 
AXS01970 
AXS01980 
AXS01990 
AXS02000 
AXS02010 
AXS02020 
AXS02030 
AXS02040 
AXS02050 
AXS0206 0 
AXS0207 0 
AXS0208 0 



C AXS0209 0 

WRITE(6,11)1,J, AMCOR ( 1, J),R1,XC(1,J),TETAC(1,J), TO, PO, VELO, AKNO, FPAXS02100 
10, DNO, P AXS021 1 0 

11 FORMAT ( ' ',2IA,5F10.3,3E12.3,F9.A,2E12.3) AXS02120 

30 CONTINUE AXS02130 

C AXS021A0 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS02150 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS02160 



non noon o non o o o o oooooo 



AXS0217 
AXS0218 
AXS0219 
AXS0220 
AXS0221 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS0222 



CALCULATE THE FLOW PARAMETERS IN THE CORE BOUNDED BY THE TWO MACH 
WAVES STARTING AT THE NOZZLE LIPS (CORNER POINTS) 

AT THE EXIT THE FLOW IS ASSUMED TO BE UNIFORM 



198 



WRITE (6,198) 
FORMAT ( ' 1 ' , ' 
APM=0 . 

BPM=0. 

CPM=0 . 



CORE 



'/) 



= 2, N1 



AMIR = 
ALFAR = 



XC(I,J) = 



AXS0223 
AXS0224 
AXS0225 
AXS0226 
AXS0227 
AXS0228 
AXS0229 
AXS0230 
AXS0231 
AXS0232 
AXS0233 
AXS0234 
AXS0235 
AXS0236 
AXS0237 
AXS0238 

(RC(I-1,J)-RC(I-1,J+1)+XC(I-1,J)*DTAN(ALFAL) +XC( 1-1 ,J+1)AXS0239 



DO 199 I 
WRITE (6,10) 

DO 199 J = 1 , N1 
IF ((I+J-l) .GT.N1) GO TO 199 
AMIL = DARSIN( 1 ./ AMCOR ( I — 1 , J ) ) 
ALFAL = PI-(TETAC(I-1, J)+AMIL) 



DARSINd . /AMCOR (1-1, J+l)) 
TETAC(I-1,J+1)-AMIR 



1*DTAN( ALFAR) )/(DTAN(ALFAL)+DTAN( ALFAR)) 

CENTR = DABS(XC( I , J ) ) 

IF (CENTR. LT.O .001) XC(I,J) = 0. 

RC( I , J ) = RC(I-1, J+1)+(XC(I, J)-XC(I-1,J+1))*DTAN(ALFAR) 



DKSI = 
DETA = 



DSQRT( (XCd, J)-XC(I-1, J + l) )XX2+(RC(I, J)-RC(I-1, J + l) )xx2) 
DSQRT( (XCd, J)-XC(I-1, J))**2 + (RC(I, J)-RC(I-1, J))**2) 



CALCULATE NOW THE PRANDTL MEYER FUNCTION AND FLOW ANGLE IN CORE. 
APM = PMC(I-1, J)+PMC(I-1, J+1)+TETAC(I-1, J+l )-TETAC( 1-1 , J) 

IF (KD.EQ.2) GO TO 151 

BPM = DSIN(AMIR)XDSIN(TETAC(I-1, J+l ) )/RC( 1-1 , J+1)*DKSI 
CPM = DSIN(AMIL)*DSIN(TETAC(I-1, J))/RC(I-1, J)xDETA 
PMC( I , J ) = (APM+BPM+CPM)/2.0 



151 



APM = PMC(I-1, J+1)-PMC(I-1, J)+TETAC(I-1, J+l )+TETAC( 1-1 , J) 
TETAC( I , J ) = (APM+BPM-CPM)/2.0 



DEFINE MACH NUMBER (BY ITERATIONS) 
INITIAL GUESS AMCOR ( I , J )= AMCOR ( 1-1 , J+l ) 



154 



160 



156 



160 



AMG = AMCOR(I-l, J+l) 

KZ = 0 

IF (KZ.GE.100) GO TO 
KZ = KZ+1 
C = AMG*AMG*A3+1 . 

D = DSQRT( AMGXAMG-1 . ) 

PMCAL = A2*DATAN(B2*D)-DATAN(D) 
DELNI = PMCAL - PMC(I,J) 

DEL = DABS(DELNI) 

IF (DEL. LT. .000002) GO TO 160 
IF (DELNI. LT. 0. ) GO TO 156 
AMG = AMG* . 999 
GO TO 154 

AMG = AMG*(1 .-DELNIXC/D) 

GO TO 154 
AMCOR( I , J ) = AMG 



CALCULATE FLOW PARAMETERS 
IF (J.GT.l) GO TO 197 
C = AMC0R(I,J)XAMC0R(I,J)*A3+1 
PRES = PSTG/ ( CXXB1 ) 

TEMP = TSTG/C 

DN = PRES/( BOLTZXTEMP) 

FP = . 7 07/( DNXCXS) 

SCALE = DKSI X DSIN(ALFAL) 

AKN = FP/SCALE 



AXS0240 
AXS0241 
AXS0242 
AXS0243 
AXS0244 
AXS0245 
AXS0246 
AXS0247 
AXS0248 
AXS0249 
AXS0250 
AXS0251 
AXS0252 
AXS0253 
AXS0254 
AXS0255 
AXS0256 
AXS0257 
AXS0258 
AXS0259 
AXS0260 
AXS0261 
AXS0262 
AXS0263 
AXS0264 
AXS0265 
AXS0266 
AXS0267 
AXS0268 
AXS0269 
AXS0270 
AXS0271 
AXS0272 
AXS0273' 
AXS0274I 
AXS0275I 
AXS027 6 1 
AXS0277 I 
. AXS027 8I 



AXS0279I 
AXS0280I 
AXS0281I 
AXS0282I 
AXS0283I 
AXS0284I 
AXS0285I 
AX50286 ( 
AXS0287 ( 
AXS0288I 



50 



c 

c 



SOUND = DSQRT(GAMAxRjxTEMP) 

VEL = SOUNDXAMCORCI, J) 

BREAKDOWN PARAMETER AS DEFINED BY 'BIRD'. 
DENSFCI,J) = PRES/(RJ*TEMP) 

DDENS = DENSF(I-1,J) - DENSF(I,J) 

COLF = 4.XCXS*DN*DSQRT(B0LTZ*TEMP/(PIXGMM)) 
P = VEL*DDENS/(SCALE*DENSFCI, J)*COLF) 

197 CONTINUE 



C 

C 

C 

C 



AXS0289 0 
AXS0290 0 
AXS02910 
AXS02920 
AXS02930 
AXS02940 
AXS02950 
AXS0296 0 
AXS0297 0 
AXS02980 
AXS02990 
AXS030 00 
AXS030 1 0 



TIME1 = SCALE/VEL 

TIME2 = 3./(4.*CXS*DN*DSQRT(B0LTZ*TEMP/(PI*GMM))) 

P = TIME2/TIME1 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS03020 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS03030 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS03040 
WRITEC6, 11)1, J,AMCORCI, J),RC(I, J),XC(I, J),TETAC(I, J ), TEMP, PRES, VELAXS03050 
1,AKN,FP,DN,P AXS0306 0 

199 CONTINUE AXSD3070 

C AXS03080 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxaXS03090 

exxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXSOSI 00 

C MATCH CORE AND FAN POINTS AXS03110 

DO 200 1=1, N1 AXS03120 

X(I,1) =XC(I,1) AXS03130 

R( I , 1 ) = RC( I , 1 ) AXS03140 

AM( I , 1 ) = AMCORC 1,1) AXS03150 

PMC 1,1) = PMCC 1,1) AXS03160 

200 TETAC 1,1) = TETACC I , 1 ) AXS03170 

C AXS03180 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS03190 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS03200 
C AXS0321 0 

C CALCULATE FLOW PARAMETERS IN REGION 2 (SIMPLE PRANDTL MEYER FAN). AXS03220 

C AXS03230 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS 03240 



GO TO 299 



C 

C 



WRITE (6,298) 
298 FORMAT ( '1', ' 

KFIN = 51 



DO 299 I = 2,N1 
WRITEC6, 10) 

DO 299 J = 2, N2 
IF (J.GT.KFIN) 

KZ = 0 

AMIL =DARSIN(1./AMCI-1, J)) 

AMIR =DARSIN(1 ./AMCI, J-l) ) 

ALFAL = PI-CTETACI-1,J)+AMIL) 

ALFAR = TETA(I, J-D-AMIR 

CHECK ANGLES AND CALCULATE COORDINATES 
ANGLE1 = PI/2. -.000001 
ANGLE2 = PI/2. +.000001 
IF (ALFAL.LE.ANGLE1 .OR. ALFAL .GE.ANGLE2) 
XCI,J) = X(I-1,J) 

GO TO 207 

201 IF (ALFAR. LE.ANGLE1 .OR. ALFAR. GE.ANGLE2) 
XCI,J) = X(I,J-1) 

GO TO 207 



205 XCI,J) = C RC I - 



AXS03250 

REGION 2V) AXS0326 0 

AXS0327 0 
AXS03280 
AXS0329 0 
AXS03300 
AXS0331 0 
AXS03320 
AXS03330 
AXS03340 
AXS03350 
AXS0336 0 
AXS0337 0 
AXS0338 0 
AXS03390 
AXS03400 
AXS0341 0 
AXS03420 

GO TO 201 AXS03430 

AXS03440 
AXS03450 

GO TO 205 AXS0346 0 

AXS0347 0 
AXS03480 
AXS03490 

1, J)-R(I, J-1)+X(I, J-1)XDTAN(ALFAR)+XCI-1, J )XDTANC ALFAXS03500 



1AL) )/(DTAN( ALFAL )+DTAN( ALFAR)) 

207 IF (ALFAL. LE. 0.000001. OR. ALFAL. GE.O 
R( I , J ) = R( 1-1 , J ) 

GO TO 213 

209 IF (ALFAR. LE. 0.000001. OR. ALFAR. GE.O 
R( I , J ) = R( I , J-l ) 

GO TO 213 



000001) GO TO 209 



000001) GO TO 211 



211 RC I , J ) = (X(I, J)-XCI,J-1))XDTAN(ALFAR)+RCI, J-l) 



AXS0351 0 
AXS03520 
AXS03530 
AXS03540 
AXS03550 
AXS0356 0 
AXS0557 0 
AXS03580 
AXS03590 
AXS03600 



51 



oo oo o on o o o o ooo o oooo 



213 DKSI = DSQRT((R(I,J)-R(I,J-1) )xx2+(X(I,J)-X(I,J-l) )xx2 ) 
DETA = DSQRT(CR(I,J)-R(I-l,J))xx2+(X(I,J)-XCI-l,J))xx2) 



IF (RCI,J).GT.O..AND.X(I,J).GT.O.) GO TO 219 
KFIN = J-l 
WRITEC6 , 12) 

12 FORMAT ( ' FURTHER POINTS ON THE CHARACTERISTICS ARE DIVERGENT’ 

219 CONTINUE 

LOCATION OF THE NEW MESH POINT HAS BEEN FOUND 



CALCULATE NOW PRANDTL MEYER FUNCTION AND FLOW DIRECTION FOR NEW POINT 
APM = PMCI, J-1)+PM(I-1, J)-TETACI-1,J)+TETACI,J-1) 

IF (KD.EQ.2) GO TO 251 

BPM = DSIN(AMIR)*DSINCTETACI, J-1))/R(I, J-1)*DKSI 
CPM = DSIN(AMIL)XDSIN(TETA(I-1, J))/R(I~1, J)xDETA 
251 PMC I , J ) = . 5# ( APM+BPM+CPM) 



APM = PMCI,J-1)-PMCI-1,J)+TETACI,J-1)+TETACI-1, J) 
TETA(I,J) = . 5X( APM+BPM-CPM) 



CALCULATE NOW MACH NUMBER FOR EACH POINT 
INITIAL GUESS AM(I,J) = AM(I-1,J) 

AMG = AMCI-1,J) 

KZ = 0 

254 IF (AMG. GT. 200 . ) GO TO 257 
D = DSQRT(AMG*AMG-1 • ) 

GO TO 258 

257 D = AMG 

258 IF CKZ.GE.100) GO TO 260 
KZ = KZ + 1 

PMCAL = A2XDATAN(B2*D)-DATANCD) 

DELNI = PMCAL - PM(I,J) 

DEL = DABSCDELNI) 

IF (DEL. LT. .000002) GO TO 260 
IF (DELNI. LT.O.)GO TO 256 
AMG = AMG* . 999 
GO TO 254 



256 IF (AMG. LT. 2000. )GO TO 2560 
KFIN = J-l 
GO TO 299 



AXS0361C 
AXS0362C 
AXS0363( 
AXS0364C 
AXS0365( 
AXS0366( 
)AXS0367( 
AXS0368( 
AXS0369( 
AXS0370I 
AXS037K 
AXS0372I 
AXS037 3( 
AXS0374I 
AXS037 5( 
AXS0376I 
AXS0377I 
AXSJ3378I 
AXS0379I 
AXS0330I 
AXS038 1 I 
AXS0382I 
AXS0383I 
AXS0384I 
AXS0385I 
AXS0386 I 
AXS0387I 
AXS0388I 
AXS0389I 
AXS0390I 
AXS0391I 
AXS0392I 
AXS0393I 
AXS0394I 
AXS0395I 
AXS0396 1 
AXS0397 I 
AXS0398I 
AXS0599I 
AXS0400I 
AXS0401I 
AXS0402I 



2560 D1 = (A3XAMGXAMG+1 . ) 

AMG = AMG*( l.-DELNIxDl/D) 

GO TO 254 

260 AM( I , J ) = AMG 

CALCULATE NOW LOCAL TEMPERATURE, PRESSURE, VELOCITY, KNUDSEN NO. 

C = AM( I , J )XAM( I , J ) XA3+1 
PRES = PSTG/(C**B1) 

TEMP = TSTG/C 

DN = PRES/CBOLTZXTEMP) 

DENSF( I , J ) = PRES/(RJ*TEMP) 

DDENS = DENSF(I-1, J-l )-DENSF( I , J) 

SCALE = DSQRT((X(I,J)-X(I-l,J-l))xx2+(R(I,J)-R(I-l,J-l))xx2) 
*x****xxxx*#x*#xxx**x***x#***xx*xxxx***x*x* 

271 FP = . 7 07/ ( DNxCXS ) 

AKN = FP/SCALE 

SOUND =DSQRT(GAMAXRJXTEMP) 

VEL = AM( I , J )xSOUND 



COLF = 4.XCXSXDNXDSQRT(B0LTZXTEMP/(PIXGMM)) 

P = VELxDDENS/(SCALExDENSF(I,J)xCOLF) 

PRINT RESULTS FOR MESH POINTS 
KL = (-l)xxl 
IF (KL.LT.O) GO TO 299 

WRITE(6,11)I,J,AM(I,J),R(I,J),X(I,J),TETA(I,J),TEMP,PRES,VEL,AKN, 



AXS0404I 
AXS0405I 
AXS0406I 
AXS0407I 
AXS0408I 
AXS0409I 
AXS0410I 
AXS0411I 
AXS0412I 
AXS0413I 
AXS0414I 
AXS0415I 
AXS0416I 
AXS0417 ( 
AXS04181 
AXS041 9( 
AXS0420( 
AXS0421 ( 
AXS0422C 
AXS0423C 
AXS0424C 
AXS0425C 
AXS0426 C 
AXS0427 C 
AXS0428C 
AXS0429C 
AXS0430C 
AXS0431C 
AXS0432C 



1FP,DN,P 
299 CONTINUE 



AXS04330 
AXS04340 

C AXS04350 
CXXX**X**XXXX*XXXXXX*XX#X*X*XXXXXXXXX*X*XX*X*XX*XX**XXX**X*X**XXXXXXXXXXAXS0436 0 
C*XXXXX*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS04370 
C AXS04380 
C MATCH ’REGION 2' AND ’REGION 3’ POINTS AXS04390 
C AXS04400 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX AXS044 10 



KFIN- 1 

= 1 , L 



L 

DO 300 J 
XXC1, J) 

RXC 1 , J ) 

AMX( 1 , J ) = 

PMX(1,J) 

300 TET AXC 1 , J ) = 



XC20, J) 
RC20, J) 
AMC20, J) 
PMC 20 , J ) 
TETAC20, J) 



AXS04420 
AXS04430 
AXS04440 
AXS04450 
AXS0446 0 
AXS0447 0 
AXS04480 



C*XXXX*X*XXXXXXX*XXXXXXXXXXXXXXXXXXXX*XXXXXXX*XX*XXXXXXXXXXXXXXXXXXXXXXXAXS04490 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS04500 
C AXS04510 
C CACULATE FLOW PARAMETERS FOR REGION 3 AXS04520 
C AXS04530 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS04540 



AXS04550 

397 FORMAT Cl',’ REGION 3V) AXS04560 

AXS0457 0 
AXS04580 
AXS04590 
AXS046 00 
AXS04610 

320 KZ = 0 AXS04620 

AXS04630 
AXS0464 0 
AXS046 50 
AXS0466 0 
AXS0467 0 
AXS04680 
AXS046 90 
AXS047 00 
AXS04710 
AXS04720 
AXS04730 
AXS047 4 0 
AXS047 50 
AXS047 6 0 
AXS0477 0 

301 AMIR = DARSINCl./AMXCI, J-1)) AXS04780 

AXS04790 
AXS04S00 
AXS04810 
AXS04820 
AXS04830 

302 IFCALFAR. LE. ANGLE 1 . OR. ALFAR.GE. ANGLE2) GO TO 305 AXS04840 

AXS04850 
AXS0486 0 

305 XX(I,J) = (RXCI-1, J)-RX(I, J-1)+XX(I, J-l)xDTAN(ALFAR)+XXCI-l, J)XDTAAXS04870 



WRITE (6,397) 

FORMAT Cl',' REGION 3V) 

DO 399 I = 2, L 
DO 399 J = I , L 
IF (J.GT.KFIN) GO TO 399 
IF CJ.GT.DGO TO 320 
WRITE (6,10) 

KZ = 0 

AMIL = DARSINC 1 ./AMX(I-1,J) ) 

ALFAL = PI-CTETAXCI-1, J)+AMIL) 

IF (J.GT.I) GO TO 301 

ALFAR = ALFAL 

TETAXC I , J ) = PIX.5 

TETAXC I , J-1 ) = PI-TETAX(I-1, J) 

XX(I,J) = 0. 

RX( I , J ) = RX( 1-1 , J )+XX( I-l,J)xOT AN (ALFAL) 

RX( I , J-1 ) = RX(I-1,J) 

XX(I,J-1) = -XX(I-1,J) 

DKSI = (RX(I,J)-RX(I-1, J))/DSIN(ALFAL) 

DETA = DKSI 

P MX (I, J-1) = PMX(I-1,J) 

GO TO 316 

AMIR = DARSINCl./AMXCI, J-1)) 

ALFAR = TETAXCI, J-D-AMIR 

IF (ALFAL. LE.ANGLE1. OR. ALFAL. GE.ANGLE2) GO TO 302 
XX(I,J) = XXCI-1,J) 

GO TO 307 

IFCALFAR. LE.ANGLE1. OR. ALFAR. GE.ANGLE2) GO TO 305 
XX(I,J) = XX(I,J-1) 

GO TO 307 



INC ALFAL ) )/CDTAN( ALFAL )+DTAN( ALFAR ) ) 



307 



309 



311 

315 



316 



IF (ALFAL. 
RXC I , J ) = 
GO TO 315 
IF (ALFAR. 
RX(I,J) = 
GO TO 315 
RXC I , J ) = 



LE. 0.000001 
RX( 1-1 , J ) 

LE. 0.000001 
RXCI, J-1) 



OR. ALFAL. GE. 0.000001) GO TO 309 



OR. ALFAR.GE. 0.000001) GO TO 311 



RXCI, J-1)+(XX(I, J)-XX(I, J-1 ))*DTANC ALFAR) 



DKSI = DSQRTC(RX(I,J)-RX(I,J-l))x*2+CXXCI,J)-XXCI,J-l))xx2) 
DETA = DSQRT((RX(I,J)-RX(I-1 ,J))xx2+CXXCI,J)-XXCI-1,J))xx2) 
CONTINUE 



IF ( RXC I , J ) . GE . 0 , 
KFIN = J-1 
WRITE (6,398) 



. AND . XXC I , J ) . GE . 0 . )G0 TO 319 



AXS04380 
AXS04890 
AXS04900 
AXS04910 
AXS04920 
AXS04930 
AXS04940 
AXS04950 
AXS0496 0 
AXS04970 
AXS04980 
AXS04990 
AXS05000 
AXS0501 0 
AXS05020 
AXS05030 
AXS05040 



351 



354 



357 

358 



398 FORMAT (’ ’ , ' FURTHER POINTS ON THE CHARACTERISTICS ARE DIVERGENT’ 
C 

C LOCATION OF THE NEW POINT HAS BEEN FOUND 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C 

C CALCULATE NOW P.M. ANGLE AND FLOW DIRECTION 

319 APM = PMX(I, J-1)+PMX(I-1,J)-TETAX(I-1,J)+TETAX(I, J-l) 

IF(KD.EQ.2) GO TO 351 

BPM = DSIN(AMIR)XDSIN(TETAX(I, J-1))/RX(I, J-1)XDKSI 
CPM = DSIN(AMIL)xDSIN(TETAX(I-l,J))/RX(I-l,J)xDETA 
PMX(I,J) = . 5x( APM+BPM+CPM) 

APM = PMX(I,J-1)-PMX(I-1,J)+TETAX(I,J-1)+TETAX(I-1,J) 

TETAXC I , J ) = .5X( APM+BPM+CPM) 

C 

C CALCULATE NOW THE MACH NUMBER FOR EACH POINT 
C INITIAL GUESS AMX(I,J) = AMX(I-1,J) 

AMG = AMX(I-1,J) 

KZ = 0 
KH = 0 
KL =0 

IF (AMG. GT. 2000.) GO TO 357 
D = DSQRT(AMG*AMG-1 . ) 

GO TO 358 
D = AMG 

IF ( KZ . GE. 50 )GO TO 360 
KZ = KZ +1 

PMCAL = A2*DATAN(B2XD)-DATANCD) 

DELNI = PMCAL-PMXd, J) 

DEL = DABSC DELNI ) 

IF (DEL. LT. .000002) GO TO 360 
IF (DELNI.LT. 0.)GO TO 356 
AMG = AMGX.98 
GO TO 354 

IF (AMG. LT. 5000.) GO TO 3560 
KFIN = J-l 
GO TO 399 
D1 = A3XAMGXAMG+1 . 

C X XX X X X XX X X X X X X X X X** X X X X X * X X X X X 
DDELNI = DELNI 

IF (AMG.GT.20.) DDELNI = DELNIX( .95XXKZ) 

CXXXXX*X**XX**XXXXX*XX*******XX 
AMG = AMGX(1 .-DDELNIXD1/D) 

GO TO 354 

360 AMX ( I , J ) = AMG 
C 

C CALCULATE NOW THE LOCAL TEMPERATURE, PRESSURE, VELOCITY, KNUDSEN 
C = AMX( I , J ) xAMX( I , J ) XA3+1 . 

PRES = PSTG/(CXXB1) 

TEMP = TSTG/C 

DN = PRES/( BOLTZXTEMP ) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
DENSF(I, J)=PRES/(RJXTEMP) 

DDENS=DENSF(I-1, J-l )-DENSF( I , J) 

FP = . 7 07/ ( DNXCXS ) 

SCALE=DSQRT((XX(I, J)-XX(I-1, J-l))xx2+(RX(I, J)-RX(I-1, J-l))xx2) 
AKN=FP/SCALE 



356 



3560 



C 

C 



Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



c 

c 

c 

c 

c 

c 

c 



-RX( I — 1 , J ) ) 



ALFALL = -ALFAL+PI 
AF = ALFALL-PI/2. 

AF = DABS(AF) 

IF (AF.LT. .000001) GO TO 370 
BF = l./DSQRTd. + DTAN(ALFALL)XX2) 

SCALE=BFX(RX(I, J-l )-DTAN( ALFALL )x(XX( I, J-l )-XX(I-l,J))- 
SCALE = DABS(SCALE) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

SOUND = DSQRT(GAMAXRJXTEMP) 

VEL = AMX( I , J ) XSOUND 
C 

COL F=4 .XCXSXDNXDSQRT( BOLTZXTEMP/ (PIXGMM)) 
P=VELxDDENS/(SCALExDENSF(I, J)xCOLF) 



) AXS0505 
AXS0506 
AXS0507 
AXS0508 
AXS0509 
AXS0510 
AXS0511 
AXS0512 
AXS0513 
AXS0514 
AXS0515 
AXS0516 
AXS0517 
AXS0518 
AXS0519 
AXS0520 
AXS0521 
AXS0522 
AXS0523 
AXS0524 
AXS0525 
AXS0526 
AXS0527 
AXS0528 
AXS0529 
AXS0530 
AXS0531 
AXS0532 
AXS0533 
AXS0534 
AXS0535 
AXS0536 
AXS0537 
AXS0538 
AXS0539 
AXS0540 
AXS0541 
AXS0542 
AXS0543 
AXS0544 
AXS0545 
AXS0546 
AXS0547 
AX30548 
AXS0549 
AXS0550 
AXS0551 
AXS0552 
AXS0553 
AXS0554 
AXS0555 
AXS0556 
AXS0557 
AXS0558 
AXS0559 
AXS0560 
AXS0561 
AXS0562 
AXS0563 
AXS0564 
AXS0565 
AXS0566 
AXS0567 
AXS0568 
AXS0569 
AXS0570 
AXS0571 
AXS0572 
AXS0573 
AXS0574 
AXS0575 
AXS0576 
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on 



PRINT RESULTS 

WRITEC6, 11)1, J,AMX(I, J),RX(I, J), 
1AKN,FP,DN,P 
399 CONTINUE 
STOP 
END 

SENTRY 



AXS0577 0 
AXS05780 

I,J),TETAX(I,J) , TEMP, PRES, V EL , AXS05790 

AXS05800 
AXS058 1 0 
AXS05820 
AXS05830 
AXS058A0 



A. 4 LIST OF MAIN SYMBOLS USED IN 'AXSYM' 



Parameter 

Name 


Physical Name 


Units 


Type 


Description 


PAMB 


ambient pressure 


pascals 


real 


ambient atmosphere 
pressure 


KD 


K dimensions 




integer 


KD=2 for two dimen- 
sional jet 

KD=3 for axisymmetric 
ring jet 


PI 


7T 




real 

constant 




BOLTZ 


Boltzmann constant 


Joule 

degree 


real 


1.38032 T(FZ3 
joules /decree 


AVOG 


Avo gadr o f s const ant 


molecules 

mol 


real 


6.0225 x 10^ b 
1/kmol 


RG 


Universal gas 
constant 


Joule 
mol .deg 


real 


8314.3 


AM£ 


Mach No. (E.P) 




real 


Mach number at exit 
plane 


TO 


Temperature (E.P) 


°k 


real 


Temperature at exist 
plane 


PO 


Pressure (E.P) 


pascals 


real 


Pressure at exit 
plane 


R1 


Cylinder radius 


m 


real 


Radius of the 

cy li ndr i c al ve hi c le 


XI 


0.5* nozzle width 


m 


real 


Half width of nozzle 


GAMA 






real 


averaged heat 
capacity ratio (jet) 


DIAM 


Molecule diameter 


m 


real 


averaged molecular 
diameter (jet) 


GM 


Molecular mass 


kg/kmol 


real 


averaged molecular 
weight (jet) 


RJ 


Gas constant 


joules 
kj .deg 


real 


gas constant (jet) 


CXS 




m z 


real 


collision cross 
section (hard sphere) 


GMM 


mass of a molecule 


kg 


real 


averaged mass of a 
molecule 


N1 






integer 


number of divisions 
(characteristics from 
the exist plane) 


N2 






integer 


number of character- 
istics in the Prandtl 
Meyer fan 


A1 






real 


Y - 1 

Y 


B1 






real 

J 


1_ 

1 ~ 1 
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A. 4 LIST OF MAIN SYMBOLS USED IN ' AXSYM' (CONTINUED) 



Parameter 

Name 


Physical Name 


Units 


Type 


Description 


A2 






real 


Ky + D/C y - l)] 


B2 






real 


1/A2 


A3 






real 


JLZl 

2 


C 






real 


M 2 y ~ 1 + 1 
2 


D,D1 ,D2 






real 


M z - 1 


PSTG 




pascals 


real 


stagnation pressure 


TSTG 




°k 


real 


stagnation tempera- 
ture 


DSTG 




kg/m^ 


real 


s tagna t ion dens i ty 


FSM 






real 


free stream Mach 
number 


FST 




°k 


real 


free stream tempera- 
ture 


FSD 




kg/m 3 


real 


free stream density 


AMIT 


PT 


radius 


real 


Mach angle ( tail of 
P.M. fan) 


AMIH 


PH 


radius 


real 


Mach angle (head of 
P.M. fan) 


PMT 


V7* 


radius 


real 


P.M, function (tail) 


PMH 


. HL _ _ 


radius 


real 


P.M. function (head) 


RAT , RAT2 , 
El 






real 


used to define a 
loparithmic division 
of the corner char- 
acteristics (50 lines 
in P.M. fan)(a linear 
division v/ould have 
resulted in concen- 
tration of character- 
istics at high Mach 
numbers ) 


DELM 






real 


difference of Mach 
numbers 


AM( I , J) 


n 




real 2-D 
array 


Mach number at 
location ( I, j ) (used 
for the corner) 


PM( I,J) 


V 


radius 




P.M. function at 
location (I, j ) (used 
for the corner) 


AMI 


u 


radius 


real 


Mach angle 


R(I,J) 




m 


real 2-D 
array 


radius (or ordinate) 
at (I,j) 


TETA( I, J) 


e 


radius 


real 2-D 
array 


flow direction at 

(u) 
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A. 4 LIST OF MAIN SYMBOLS USED IN ' AXSYM' (CONTINUED) 



Parameter 

Name 


Physical Name 


Units 


Type 


Description 


ALFAL 




radians 


real 


angle of a left run- 
ning characteristics 


ALFAR 




radians 


real 


angle of a right run- 
ning characteris tics 


DENSF(I.J) 




kg/m^ 


real 2-D 
array 


gas density at point 
IJ 


AMCOR( I, J) 






real 2-D 
array 


Mach number at mesh 
points in region (1) 
(region 1 is bounced 
by the two tail 
charac teris tics 


TETAC( I , J) 




radians 


real 2-D 
array 


flow direction at 
mesh points (region 
1) 


PMC( I,J) 




radians 


real 2-D 
array 


P.M. function at mesh 
points (region 1) 


XC(I,J) 




m 


real 2-D 
array 


X coordinate at mesh 
points (region 1) 


RC( I,J) 




m 


real 2-D 
array 


radius (ordinate) at 
mesh points (region 
1) 


AHIL 




radians 


real 


Mach angle far left 
running character- 
is tics 


AMIR 




radians 


real 


Mach angle far right 
running character- 
is tics 


DKSI 




m 


real 


distance between mesh 
points 


DETA 


dn 


m 


real 


distance between mesh 
points 


DELHI 


dv 


radians 


real 


P.M. differential 


AMG 




radians 




Mach number (used for 
iterative calcula- 
tion) 


PRES 








pressure 


TEMP 








temperature 


DN 








number density 


FP 








mean free path 


AKN 








Knuds en number 


SOUND 








speed of sound 


VEL 








absolute local 
velocity 


DDENS 








density difference 
between two points a- 
long a steamline 


COLF 








collis ion frequency 


P 








breakdown parameter 
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A. 5 AXSYM PROGRAM USER 1 S GUIDE 



1. Input data: 

ambient pressure PAMB 

Mach number (Exit plane) AMO 

Temperature (Exit plane) TO 

Pressure (Exit plane) PO 

Half width of nozzle XI 

Radius of nozzle ring R1 

Specific heat ratio of jet gas GAMA 

average molecular diameter (jet) DIAM 

average molecular weight GM 

2. Options for flow geometry 

two dimensional flow KD=2 

ring jet KD=3 

default condition KD=3 

3. Resolution of mesh points 



to change the resolution of the mesh points 
a - change N1 and N2 as necessary 

b - change DIMENSIONS 1 according to new values of and N2 

c - define distribution of Mach lines in the Prand tl-Meyer fan as required 
(program lines 126-131) 

4. Execution commands: 

After copying program into USER* S FILE: 

WAT F IV AXSYM * (XTYPE 

The program will run on user’s terminal under WATFIV. A soft copy of 
program listing and output listing will be stored in user’s disk named 
AXSYM LISTING. 
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5. Hard Copy 



PRINT AXSYM LISTING. 

6. Program Outputs. 

All necessary outputs are automatically listed by the program. 

Figure (20) and Figure (21) show the resulting mesh of characteristics 
calculated for an altitude of 200 km for 110=4 and 110=2 respectively. 

In these figures we also show some isotherms and the limit where the 
breakdown parameter equals 0.05. These lines are plotted (manually) using 
interpolation procedures. Data along the breakdown line is input data for 
the molecular flow. 



60 



APPENDIX B SIMUL PROGRAM 



B.l DATA ORGANIZATION 

Because of the large number of molecules, cells, regions and sectors in 
the simulation and the large number of data related to each molecule, each 
cell and region to be stored, special precautions should be taken in order not 
to overflow the available computer memory. 

The following data organization was used in SIMUL. Figure 22 shows the 
geometry of one sector. 



61 




Figure 22. Definition of sector geometry and cell volume. 

cell volume = v(I) = R(I) * DDALFA x DR x RSI * DFI. 

v ( I) = R( I) ♦ RSI(I) 
v( 1 ) R( 1 ) • RSI 

v(l) is the volume of smallest cell in a region. 
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1. Tables of Molecules and Their Parameters 



PI (L ,N1M) 


Light molecules of the jet 


P2 ( L , N2M) 


- Heavy molecules of the jet 


P3(L,N3M) 


- Ambient gas molecules (not used in the present program) 



N1M, N2M, N3M - maximum number of molecules in simulation. Number of 





active molecules may be smaller or equal to (N<I>M). 


L-1,2,3 


- cartesian components of velocity v x , Vy, v z [m/s] 


L=4 


radial coordinate [meters] 

if r=-99 the particular molecule is inactive 


L=5 


- angular coordinate [radians]. 



This table is generated each time the simulation is initiated for a region. 
That means, the same group of molecules (as stored in the computer memory) is 
used to simulate the flow in all regions in the computation domain. 
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2. C(M,I,j) Table of Cells (In a Region)-Real Data 



I = 1,10 
j = 1,10 
M = 19 
M = 20 
M = 1,9 

M = 10-18 



radial index of the cell in a region 
angular index of the cell in a region 
radial coordinate of cell center 
angular coordinate of cell center 

time parameter for collisions of different species in a 
cell 

maximum relative velocity expected for collisions of 
different species in a cell 
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3 . IP(N1A+N2A+N3A) Table of the addresses of the active molecules 

arranged in order of their species and in the order of their cells 
IC(N,I,j) - table of cells (in a region) integer data 
I = 1,10 - radial index 

j = 1,10 - angular index 

N = 1 - number of molecules (spec 1) 

N = 2 - number of molecules (spec 2) 

N = 3 - number of molecules (spec 3) 

N = 4 - (starting address - 1) of molecules as ordered in (IP) 
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4. Reg(N,kR,kS) Data table for a specific region (real) 



kR = 1,10 
kS = 1,20 
N = 1 

N = 2 
N = 3 
N = 4 
N = 5 



- index of a region in a sector 

- index of the sector 

- Dfl = differential angle (axisymmetric) 

DF1 is a weighting factor 

- DN1 = number density (species 1) 

DN2 = number density (species 1) 

- VOL 1 = actual volume of smallest cell in kR 

- AREA1 INPUT area of smallest cell in kR 
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5 . Region Geometry and Input Flux 



R( 10) 
A( 10) 

V0L(10) 
Ml (10) 

Fl( 10) ,F2(10) 



polar radius of a cell in a region 

input area of a cell 

input area of smallest cell 

volume of a cell 

volume of smallest cell 

initial number of molecules in cells 

(equal number of molecules of either jet species) 

input flux through high pressure starting line 

(spec 1), (spec 2) 
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6 . Input Flux From High Pressure Starting Line 



FWPl(N,I,j) 

FWP2(N, I, j ) 
+ t'M‘ 



species 

positive (input molecules) 
wes t 
flux 



I = 1,10 



j « 1,10 

N = 1 
N = 2 
N = 3 
N = 4 



number of the cell along the starting line in a region 

(j) 

number of the region along the starting line 

molecular flow for a given DF1 

mean molecular velocity Vx 

mean molecular velocity Vy 

mean input gas temperature 



7 • Output Flux 



FNNl(4,NAD,kR) 



FSNl(4,NAD,kR) 



FNN2(4,NAD,kR) 

ft 



FSN2(4,NAD,kR) 

ft 



(negative (output) (negative (output) 

I I 

north south 

kR - number of region in the sector 
NAD - angular location 

The first index include the same parameters as (FWP1) 



FEN1 ( 4 ,NRD) 
FEN2(4,NRD) 
(negative 

I 

eas t 



FWN1 (4 ,NRD) 

FWN2(4,NRD) 

ft 

negative 
wes t 
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NRD ■+■ radial location 



FEN1 ,FEN2 ,FWN1 ,FWN2 are necessary for iterations within one sector. 

8 . Sampling of Output Flux from a Sector 

After averaging they are transferred to F0E1, F0E2, F0W1 , F0W2. 

F0W1 (4 ,kC ,kR,kS) 

output flow to the west 

FOW2(4,kC,kR,kS) 

F0E1 ( 4 ,kC ,kR,kS) 

output flow to the east 

FOE2(4,kC,kR,kS) 

kC - radial location (cell) of flow in a region 
kR - index of a region in a sector 
kS - index of the sector 

The first index (4 parameters) include the same parameters as FEP1( ). 
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B.2 MOLECULAR SIMULATION FOR A GIVEN REGION 



After we define the geometry of the whole sector resulting from the 
region geometry and cell geometry, we may start with the molecular simulation. 
This includes : 

- initial setting of molecules in cells 

- molecules are moved according to the time increment DTH 

- new molecules are generated according to input (or output) flows 

- collisions calculations 

- integration of flow parameters for average parameters calculations 

- repetition of the whole procedure as long as necessary to obtain 
reasonable statistical averaging 

- calculate averages and flow weighting. 

These routines are the core of the program and must be repeated for all 
regions in a sector and (or all sectors in the domain where the 
collisions are significant). 

In the following sections we bring a detailed description of this part of 
the program. 

1 • Initial Setting of Molecules in Cells 

a. The initial number of available molecules in a simulation is 
larger than the number of active molecules (PI (data, number of 
molecules), P2(data, number of molecules) are the vectors used for 
species 1 and 2) 

b. an inactive molecule is defined as 

[PI (4 ,N) or P2(4 ,N) ] = -99 
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c. calculation of number of molecules to be set in each cell 

d. calculation of cell coordinates 

e. deactivation of all available molecules in simulation 

f. definition of molecules coordinates 
PKAjN), P2(4,M) are polar radiuses 

P1(5,N), P2(5,M) are angular coordinate in radians 

All molecules in a cell are set at random locations within the 
cell . 

g. definition of molecular velocities 
P1(1,N), P2(1,M) velocity in X direction 
P1(2,N), P2(2,M) velocity in Y direction 
P1(3,N), P2(3,M) velocity in Z direction 

Thermal velocities are random function of temperatures and are 
added to the mean velocity as defined at initial boundary ALFA Q of 
the region. 

As the thermal velocity has a Boltzmann distribution the thermal velocity 
setting is based on rejection-acceptance methods (for more details see Bird 
[4] Appendix D) . 

h. reset collision timers and relative velocity 
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i. reset general time counter: Time = 0 

3 . The Simulation 

a. Move all molecules according to their velocity (V x , Vy) and find 
their new coordinates 

Note - A routine designed to calculate the collisions with the 
wall was included in the program; if the region (or sector) is 
bounded by the solid wall, the collision is calculated - resulting 
new velocities and directions and counted for wall flux 
calculations. If the program is stopped at an angle where the 
flow becomes collisionless, this routine becomes irrelevant and 
other type of calculations should be designed. 

b. Output flow counting: all molecules that leave the region are 

counted and stored in specific vectors which are used as 

inputs to other regions. The output vectors are (X represent 1 or 
2 for the two species in the program (FSNX)) 

FSNX(l,j,kR) -► "south" boundary 
FNNX(l,j,kR) > "north" boundary 
FENX(1,I) -> "east" boundary 

FWNX(1,I) * "west" boundary 

I represents the radial location index of the cell 
j represents the angular location index 
kR is the index of the region within a sector 
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Note - all molecules that move to 



( j =NAEH-1 ) ( I=NRD+1 ) are placed in FENX(1,NRD) 

(j = -1)(X=NRLH-1) are placed in FWNX(1,NRD) 

(j=NAEH-l)(I= -1) are placed in FSNX(l,l,kR) 

(j= -1 ) ( 1= -1) are placed in FWNX(l,l,kR) 

This was done only for simplification reasons, 
c. Generation of new molecules due to input flows. Through the four 
sides E,N,W,S, of a region, molecules are allowed to enter the 
region according with the flows coming from the neighbouring 
regions : 

for "W M side of the region in sector 1 

FWPX(P,I,kR) 

for other sectors 

FOEX(P,I,kR,kS-l) 

for "E M side of any region 

FOWX(P,I,kR,kS+l) 

for "N" side of any region 

FSNX( P , j ,kR+l) 

for M S" side of the region 

FNNX( P , j ,kR-l ) 

The first parameter of all these arrays represent: 

P = 1 - number flux (real number) 

P = 2 - velocity component - (V x ) 

P = 3 - velocity component - (Vy) 

P = 4 - gas temperature 
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Note 1 - because every region has a different size of angle DFI, all 
fluxes have to be adjusted accordingly. 

Note 2 - input fluxes are calculated, adjusted and stored as real 
numbers. The number of input molecules are by definition integers. 

In order not to ’’loose" molecules, the number of input molecules is 
increased by 1 on a random basis. (The average of many runs will 
result in the accurate average input flow.) 

New molecules are set at random locations on the boundary of the specific 
cells and at random time within DTM. Then each molecule is allowed to enter 
the region according to its initial coordinate and velocity. At the end of 
the time interval the new location and velocity is stored in molecule array 
PI or P2. If DTM is chosen to be too large and cell size is small (total 
region size too small) some molecules may cross the region and will not be 
counted in the simulation of the specific region. In order not to "loose” 
molecules : 

(a) DTM should be decreased 

(b) count these molecules as output fluxes from the specific 
region. (This is recommended only if there is no other 
choice.) 

Note - Because the arrays of input flows store only averaged data for 
the molecules, the thermal velocity of each new molecule is calculated 
according to the Boltzmann distribution as a function of the averaged 
temperature . 

d. Rearrangement of molecules in cells. Before collisions are 

calculated all simulated molecules which have been let to move and 
generated have to be rearranged and recounted for each cell. 
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The array IC(k,I,j) contains integer data for each cell 

(I,j) 

k a 1 - is the number of molecules of species 1 

k = 2 - is the number of molecules of species 2 

k = 3 - is the number of molecules of species 3 (not used) 

k - 4 - is the (address-1) of the first molecule in the cell 

related to the vector IP(M) 

Vector IP(M) contains the list of the simulated molecules arranged 
in the order of species in cells and cells respectively* The 
following is a graphic description of IP(M) 







Figure 23. Vector IP(M) 
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B.3 SIMUL Program Flowchart 

A simplified flowchart for the Monte Carlo simulation of the molecular 
flow is shown in Figure 24. The program is designed to solve the ring 
axisymmetric jet flow, however, minor changes may be done to enable a 
different geometry. 
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Figure 24. 
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B.4 SIMUL Program Listing 



C PROGRAM SIMUL SIM00010 

C THIS PROGRAM IS DESIGNED TO CALCULATE THE MOLECULAR FLOW OUTSIDE THE SIM00020 
C CONTINUUM REGION FLOW OF A HIGHLY UNDEREXPANDED AXISYMMETRIC RING JET.SIM00030 
C RESULTS FOR THE CONTINUUM FLOW MAY BE OBTAINED FROM 'AXSYM' PROGRAM SIM00040 
C WHICH GIVES THE 'METHOD OF CHARACTERISTICS' ISENTROPIC SOLUTION. 

C THE BOUNDARY BETWEEN THE CONTINUUM AND MOLECULAR FLOW IS DEFINED BY 
C THE BREAKDOWN PARAMETER 'P' AS PROPOSED BY G. A. BIRD. 

C THE 'MOLECULAR DOMAIN' IS DIVIDED INTO POLAR DIVISIONS MAKING A SET 
C OF SECTORS. EACH SECTOR IS SUDIVIDED INTO 10 REGIONS 
C NO APRIORY INFORMATION ABOUT THE GEOMETRY OF THE DIFFERENT 
C REGIONS IS AVAILABLE, THEREFORE, MANUAL INTERVENTION MAY BE 
C REQUIRED WHEN MOVING FROM ONE SECTOR TO THE OTHER. AN ACCEPTABLE 
C GEOMETRY WILL RESULT A REASONABLE NUMBER OF SIMULATED MOLECULES. 

C EACH REGION IS SUBDIVIDED INTO A NUMBER OF CELLS WITH A GEOMETRY 
C DEFINED BY A POLAR MESH. A NUMBER OF MOLECULES IS SET IN EACH CELL 
C PROPORTIONAL TO THE CELL VOLUME. THE BOUNDARY CONDITIONS FOR EACH 
C REGION REQUIRE INPUT AND OUTPUT FLUX OF MOLECULES. NO APRIORY 
C INFORMATION ON THE FLUX IS AVAILABLE. IT WILL BE CALCULATED IN 
C AN ITERATIVE MODE. 

C IF THE BOUNDARY CONDITIONS FOR ALL CELLS ARE CONSTANT THE NUMBER 
C FLUX INTO EACH CELL IS PROPORTIONAL TO CELL WALL AREA. 

C TO DECREASE THE ERROR WHEN INTRODUCING MOLECULES -(INTEGER NUMBER)- 
C DUE TO THE INPUT FLUX -(REAL NUMBER), AN ADDITIONAL MOLECULE IS 
C GENERATED ON THE BASIS OF RANDOM NUMBERS SUCH THAT THE AVERAGE OF A 
C LARGE NUMBER OF SAMPLINGS WILL EQUAL THE REAL INPUT FLUX. 



SIM00050 
SIM00060 
SIM00070 
SIMQ0080 
SIM00090 
SIMOOIOO 
SIM00110 
SIM00120 
SIM00130 
SIM0D140 
SIM00150 
SIM00160 
SIM00170 
SIM00180 
SIM00190 
SIM00200 
SIM0021 0 
SIM00220 
SIM00230 
SIM00240 
SIM00250 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX SI MO 0 26 0 



THE MONTE CARLO SIMULATION 
THE JET GAS IS COMPOSED OF TWO SPECIES OF MOLECULES 
AMBIENT GAS IS REGARDED AS ONE SPECIES 
THE MOLECULAR MODEL IS - 'HARD SPHERE MOLECULE' 
NETWORK DEFINITION 



SIM00270 

SIM00280 

SIM00290 

SIM00300 

SIM00310 



THE MAXIMUM RADIUS (POLAR) OF THE DOMAIN IS ASSUMED TO BE RP=15 M. SIM00320 
THE ANGLE OF THE BOUNDARY BETWEEN CONTINUUM AND MOLECULAR FLOW SIM00330 

AND THE SOLID WALL IS SIM00340 

ALFA (CALCULATED IN PROGRAM ' AXSYM’ ) SIM00350 

DEFINE A POLAR SECTOR WITH A RADIUS 'RP' AND AN ANGLE OF SIM00360 

DALFA =5*MFP*NAD/RP SIM00370 

THIS SECTOR IS SUBDIVIDED INTO A NUMBER OF RADIAL DIVISIONS SIM00380 

MAKING 'N' REGIONS FOR SIMULATION CALCULATIONS FOR EACH 'DALFA'. SIM00390 
EACH REGION IS DEVIDED INTO SIM00400 

NAD -ANGULAR DIVISIONS (15) WITH AN ANGLE OF DDALFA=5*MFP/RP SIM00410 

NRD -RADIAL DIVISIONS (10) SIM00420 

MAKING NADXNRD CELLS SIM00430 

THE SMALLEST CELL CONTAINS 'MIN' MOLECULES SIM00440 

MIN = 15 SIM00450 

TOTAL NUMBER OF ACTIVE MOLECULES IN A REGION IS LIMITED TO 6000 SIM00460 

(3000 MOLECULES OF EACH SPECIES.) SIM00470 

THE WIDTH OF A CELL DEFINED BY ANGLE ' DFI ' MAY NOW BE EVALUATED. SIM00480 

(MIN/NUMBER DENSITY = ACTUAL CELL VOLUME) SIM00490 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM00500 
C DEFINE PARAMETERS BLOCK 1SIM00510 



COMMON IX 

DIMENSION P1(5,3000),P2(5,3000),IP(6000),C(20,10,15), 
*IC(4,10,15),REG(5,10,20) 

DIMENSION DFI ( 10 ) , DN( 1 0 ) 

DIMENSION R(10),A(10),VOL(10),M1(10),F1(10),F2(10) 

DIMENSION SS(2,10,10) 

DIMENSION FEN 1(4, 10), FWN1 (4,10),FNN1(4,15,10),FSN1(4,15,10) 
DIMENSION FEN2(4, 10), FWN2(4, 10),FNN2(4, 15, 1 0 ) , FSN2( 4, 15, 10) 
DIMENSION FW»1(4,10,10), FWP2(4, 10,10) 

DIMENSION FOE1(4,10,10,20),FOW1(4,10,10,20) 

DIMENSION F0E2(4, 10,10,20) , F0W2( 4, 1 0 , 1 0 , 20 ) 

DIMENSION NM(3),VRC(3),TOC(3,3),SPEC(3,5),NCOL(10,15) 

P1,P2,P3 CONTAIN INFORMATION ON UP TO M SIMULATED MOLECULES 

P1(1,N),P1(2,N),P1(3,N) ARE U,V,W VELOCITY COMPONENTS (CARTESIAN) 
P1(4,N),P1(5,N) ARE R AND TETA COORDINATES (POLAR) 



SIM00520 
SIM00530 
SIMO 0540 
S IMO 0550 
SIMO 056 0 
SIM0057 0 
SIM00580 
SIM00590 
SIM00600 
SIM00610 
SIM00620 
SIM00630 
SIM00640 
SIM00650 
SIM00660 
SIM 00670 
SIM00680 
SIM00690 



IP ( M ) ARE THE M MOLECULES ARRANGED IN THE ORDER OF THEIR CELLS 
C( 20 , 1 » J ) CONTAINS INFORMATION ON UP TO IXJ CELLS 
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM00700 

SIM00710 

ALFA IS THE ANGLE WHERE THE BREAKDOWN PARAMETER EQUALS .05 . SIM00720 
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C FPM IS THE MEAN FREE PATH AT ALFA. SIM00730 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM00740 
C SET GENERAL CONSTANTS BLOCK 2SIM00750 



IX 

PI 

BOLTZ 

PO 

TO 

AVOG 
RG 



529814367 
3.141593 
1 .38044E-23 
101325. 

273. 

2 . 68699E+25 
8314. 



SIM00760 
SIM0077 0 
SIM00780 
SIM00790 
SIM00800 
SIM00810 
SIM00820 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM00830 
C SET PROGRAM CONSTANTS BLOCK 3SIM00840 



I SPEC=2 
R1 = 2.5 
DR= . 15 
RP=15. 

NUMBER OF SIMULATED 
NM0L1=3000 
NM0L2=3000 
NM0L3=0 
MIN = 15 



MOLECULES IN SMALLEST CELL (DIFFERENT SPECIES) 



15 

15 

DIVISIONS 

15 



IN A SIMULATED REGION 



C 

C 

C 

c 



N1M = 

N2M = 

NUMBER OF 
NAD = 

NRDS=1 0 
NRD=1 0 
NRD1= 3 
NIS =2 

R1 IS THE RADIUS OF THE CYLINDER (WALL) 

R( I ) IS THE RADIAL COORDINATE OF CELL(I) 

VOL ( I ) IS THE RATIO BETWEEN VOLUMES OF CELL(I) AND CELL(l) 

A( I ) IS THE RATIO BETWEEN INPUT FLOW AREAS OF CELL(I) AND CELL(l) 

DR IS THE RADIAL SIZE OF A CELL (RADIAL INCREMENT) 

INPUT THE FOLLOWING AS 'DATA 1 OR 'READ' STATEMENTS xxxxxxxxxxxxxxxx 
TETA =1.2 

TETA IS THE FLOW DIRECTION (RADIANS) ON THE BREAKDOWN LINE AS FOUND 
FROM THE AXSYM PROGRAM. VO IS THE FLOW VELOCITY (M/SEC) 

VO = 2100 

VOX = VOXCOS(TETA) 

VOY = VOXSIN(TETA) 

TWALL=300 . 

DNO 1 =1.1 E21 
DN02=1 . 1E21 
DTM=1 . 5E-6 

SET DATA FOR THE DIFFERENT SPECIES 
MOLECULAR MASS 

SPEC(1,1)=4./AV0G 
SPEC(2,1)=40./AV0G 
SPEC (3,1) =29 ./AVOG 
MOLECULAR DIAMETER 
SPEC(1,2)=2. 19E-10 
00E-10 



SPEC(2,2)=4 , 
SPEC( 3,2)= 
SPEC( 1,3)= 
SPEC( 2, 3 ) = 
SPEC( 3,3)= 



ADDITIONAL DATA IF REQUIRED 



SIM00850 
SIM0086 0 
SIM0087 0 
SIM00880 
SIM00890 
SIM00900 
SIM00910 
SIM00920 
SIM00930 
SIMO 0940 
SIM00950 
SIM00960 
SIM0097 0 
SIM00980 
SIM00990 
SIMO 1000 
SIM01010 
SIM01020 
SIM01030 
SIM01040 
SIM01050 
SIM01060 
S I MO 1 07 0 
SIMO 1 08 0 
SIM01090 
SIMO 1100 
SIMO 1 1 1 0 
SIM01120 
SIM01130 
SIM01140 
SIM01150 
SIMO 1160 
SIMO 1 17 0 
SIM 01180 
SIM01190 
SIM01200 
SIM01210 
SIM01220 
SIM01230 
SIM01240 
SIM01250 
SIM0126 0 
SIM01270 
SIM01280 
SIM01290 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM01300 



C THERMAL VELOCITIES AT WALL TEMPERATURE 
VWM1 =SQRT ( 2 . XB0LTZXTWALL/SPEC( 1 , 1 ) ) 
VWM2=SQRT(2.XB0LTZXTWALL/SPEC(2,1) ) 
C 



SIM01310 

SIM01320 

SIM01330 

SIM01340 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM01350 



INITIALISATION. .RESET ALL SAMPLING VARIABLES 
DO 80 1=1, NRD 
DO 80 JR=1 , 1 0 
DO 80 JS=1 , 20 
DO 80 KPAR=1 , 4 
F0E1 ( KPAR, I,JR,JS)=0. 

F0E2(KPAR, I, JR, JS)=0 . 

F0W1(KPAR,I, JR, JS)=0 . 

80 F0W2(KPAR,I, JR,JS)=0. 



SIM01360 
SIM0137 0 
SIM01380 
SIM01390 
SIM01400 
SIM01410 
SIM01420 
SIM01430 
SIM01440 
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CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0145C 
C SET INPUT PARAMETERS BLOCK 3 SIM0146C 



ITER=1 

C RETURN TO 3000 FOR ADDITIONAL ITERATIONS 
3000 KS=1 

ALFA=1 .3 
FPM= .001 
TEMP=40 . 

VTER1=SQRT(2.*B0LTZ*TEMP/SPEC(1, 1)) 
VTER2=SQRT(2.*B0LTZ*TEMP/SPEC(2, 1)) 
IF (ITER.GT.l)GO TO 2005 
DO 2001 1=1,10 
REG(2, I , I ) =DN01 
REG( 3 ,1,1 ) =DN02 
DO 2001 J=l,10 
FWP1 C 2, 1 , J ) =V0X 
FWP2( 2, 1 , J ) =V0X 
FWP1 ( 3, I , J ) =V0Y 
FWP2(3, I, J)=V0Y 
FNP1 ( 4,1, J )=TEMP 
2001 FWP2(4, I, J)=TEMP 
2005 CONTINUE 
C 



SIM0147C 

SIM0148C 

SIM0149C 

SIM0150C 

SIM0151C 

SIM0152C 

SIM0153C 

SIM0154C 

SIM0155C 

SIM0156C 

SIM0157C 

SIM0158C 

SIM0159C 

SIM0160C 

SIM016K 

SIM0162( 

SIM0163( 

SIM0164( 

SIM0165< 

SIM0166( 

SIM0167I 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0168( 



SIM0169I 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0170< 
C DEFINE SECTOR GEOMETRY BLOCK 4 SIM0171I 



C KR IS THE INDEX FOR A REGION IN ONE SECTOR 
C RETURN TO 1000 FOR NEXT SECTOR 
DAL FA = 0 . 

1000 KR=1 

C RESET SAMPLING OF FLOW VARIABLES (N AND S) 

DO 81 1=1, NAD 
DO 81 JR=1 , 1 0 
DO 81 KPAR=1 , 4 
FNM1 C KPAR, I , JR ) =0 . 

FNN2CKPAR, I , JR) =0 . 

FSN1CKPAR, I, JR)=0 . 

81 FSN2C KPAR, I , JR) =0 . 

ALFA IS THE STARTING ANGLE OF THE PRESENT SECTOR. PREVIOUS DALFA 
SUBTRACTED FROM PREVIOUS ALFA (PREVIOUS SECTOR) GIVES THE NEW ALFA. 
ALFA=ALFA-DALFA 
DALFA = 5 . *FPM* FLOAT (NAD)/RP 
DALFA IS THE ANGLE OF THE SECTOR 

I F( DAL FA. GT. ALFA/2. ) DAL FA=AL FA/2. 

I F( DAL FA . GT . ( FPM#7 5/RP) ) DALFA = AL FA 
I F( ALFA . LE . 0 . )GO TO 2000 
DDALFA = DAL FA/ FLOAT (NAD) 

DDALFA IS THE ANGLE OF A SINGLE CELL IN THE SECTOR 
ALFAJ=ALFA-DALFA/2. 



C 

C 



SIM0172I 
SIM0173I 
SIM017 4( 
SIM0175I 
SIM0176I 
SIM0177( 
SIM0178( 
SIM0179I 
SIM0180( 
SIM0181I 
SIM0182I 
SIM0183I 
SIM0184I 
SIM0185I 
SIM0186I 
SIM0187( 
SIM0188( 
SIM0189( 
SIM0190( 
SIM019K 
SIM0192( 
SIM0193( 
SIM0194( 
SIM0195( 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0196( 



C DEFINE REGION KR IN SECTOR 

C RETURN TO 2000 FOR THE NEXT REGION KR 
2000 CONTINUE 

C RESET SAMPLING OF FLOW VARIABLES PER REGION 
DO 82 1=1, NRD 
DO 82 KPAR=1 , 4 
FENKKPAR, I) = 0. 

FEN2( KPAR, I ) =0 . 

FWN1 ( KPAR , I ) =0 . 

82 FWN2( KPAR, I ) =0 . 

MT = 0 

NRD=NRDS 

IF( KR . LT . 2)NRD=NRD1 
DO 100 I = 1 , NRD 

C POLAR RADIUS MEASURED FROM THE NOZZLE LIP 
R( I ) =(FL0AT(I)-.5)*DR 
IF(KR.EQ.2)R(I)=R(I)+FL0AT(NRD1)XDR 
IF(KR.GE. 3)R(I)=R(I )+( FLOAT (NRD1+(KR-2)XNRD) )XDR 
C RSI IS THE RADIUS MEASURED FROM THE AXIS OF SYMMETRY 
RSI=R(I)+R1/SIN(ALFAJ) 



BLOCK 5 SIM0197( 
SIM0198( 
SIM0199( 
SIM0200( 
SIM020K 
SIM02021 
SIM0203C 
SIM0204( 
SIM0205C 
SIM0206 ( 
SIM0207C 
SIM0208C 
SIM0209C 
SIM0210C 
SIM0211C 
SIM0212C 
SIM0213C 
SIM0214C 
SIM0215C 
SIM0216 C 
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IF(I.EQ.1)RS1=RSI 

A( I ) = RSI/RSI 

VOL ( I ) = RCI)XRSI/CRC1)XRS1) 

Ml ( I ) = MINxVOLCI) 

MT = MT+M1CI) 

MIC) IS THE INITIAL NUMBER OF MOLECULES IN EACH CELL 
DN1=REGC2,KR,KS) 

DN2=REGC3,KR,KS) 

REG(1,KR,KS)=FL0AT(MIN)/(DN1*R(1)XRS1*DR*DDALFA) 

DFI ( REGC 1 , KR, KS) ) IS THE SPHERICAL ( AXISYMMETRIC) ANGLE 

THIS ANGLE HAS DIFFERENT VALUES FOR EACH KR, THEREFORE IT IS A 

WEIGHTING FACTOR. +=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+++++++++++++ 



C. 

C 



CFNN1,FNN2,FSN1,FSN2,FEN1,FEN2,FWN1,FWN2) X DFI 

CF0E1,F0E2,F0W1,F0W2) x DFI 

DAI = RS1XDRXREG(1,KR,KS) 

DAI IS THE INPUT AREA OF CELL(l) 

FN1 = V0XDA1XREG(2,KR,KS)XSINCALFA-TETA) 

FN2 = VOXDA1XREGC 3, KR, KS )XSIN( ALFA-TETA) 

TETA IS THE ANGLE BETWEEN FLOW DIRECTION AND THE WALL 
F1CI) = FN1XACI) 

F2CI) = FN2XACI) 

99 FWP1 ( 1,I,KR)=F1CI) 

100 FWP2C 1,I,KR)=F2CI) 



DEFINE ACTUAL VOLUME AND INPUT AREA OF SMALLEST CELL IN REGION 
REG(4,KR,KS)=R(1)XDDALFAXDRXRS1XREGC1,KR,KS) 

REGC 5, KR, KS ) =DA1 

C FN1 IS THE NUMBER FLUX TO THE SMALLEST CELL CREAL NUMBER) PER SECOND 
C INPUT NUMBER OF MOLECULES CINTEGER) WILL BE INTEGRATED TO MAKE AN 
C AVERAGE OF FN1 . 

C F1CI) IS THE INPUT FLUX TO CELLCI) 

C FOLLOWING ARE THE REGION BOUNDARIES 

102 RMIN=RC1)-.5XDR 
RMAX=RMIN+DRXFLOATCNRD) 

TMAX=AL FA 
TMIN=ALFA-DALFA 

DO 9102 I = 1 , NRD 

WRITE C6,9101)RCI),ACI),VOLCI),M1CI),F1CI) 

9101 FORMAT C ' • , 3F13 .5, 110, E1S .5) 

9102 CONTINUE 
MTT - MTXNAD 
WRITEC6 , 103 )MT ,MTT 

103 FORMAT C INITIAL NUMBER OF SIMULATED MOLECULES IS= ',15,' PER 
XCIES PER DDALFA, TOTAL NUMBER IS’, 15) 



SIM0217 0 
SIM0218 0 
SIM02190 
SIM02200 
SIM02210 
SIM02220 
SIM02230 
SIM022A0 
SIM02250 
SIM02260 
SIM0227 0 
SIM02280 
SIM02290 
SIM02300 
SIM02310 
SIM02320 
SIM02330 
SIM023A0 
SIMQ2350 
SIM02360 
SIM0237 0 
SIM0238 0 
SIM02390 
SIM02A00 
SIM02A10 
SIM02A20 
SIM02A30 
SIM02AA0 
SIM02A50 
SIM02A6 0 
SIM02A70 
SIM02A8 0 
SIM02A90 
SIM02500 
SIM0251 0 
SIM02520 
SIM02530 
SIM025A0 
SIM02550 
SIM0256 0 
SIM0257 0 
SIM0258 0 
SIM02590 
SPESIM026 00 
SIM026 1 0 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM02620 

C DEFINE INPUT FLOWS TO KR W,E,N,S BLOCK 6 SIM02630 

C SIM026A0 

C SIM02650 

C SIM0266 0 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM02670 



C CALCULATE INITIAL NUMBER OF MOLECULES IN CELLS BLOCK 7 SIM02680 

C SET INITIAL STATE OF GAS " 8 SIM02690 

C SIM02700 

DO 150 I =1 , NRD SIM027 1 0 

DO 150 J= 1 , NAD SIM02720 

C SIM027 30 

C SET SIMULATED MOLECULES IN THEIR CELLS SIM027A0 

ICC1,I,J) = Ml C I ) SIM02750 

ICC 2 , 1 , J ) = Ml C I ) SIM02760 

ICC 3 , 1 , J ) = 0 SIM0277 0 

C SIM0278 0 

C SET CELL COORDINATES SIM02790 

CC 19 , 1 , J ) = RCI) SIM02800 

CC 20 , 1 , J ) = ALFA-CFL0ATCJ)-.5)XDDALFA SIM02810 

, 150 CONTINUE SIM02820 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX SIM02830 
C SIM028A0 

C DEACTIVATE ALL MOLECULES SIM02850 

DO 170 N = 1 , NMOL 1 SIM02860 

PI C A , N ) = -99. SIM0287 0 

170 P2C A, N) = -99. SIM02880 
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no o no o o no 



c 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
SET INITIAL STATE OF THE GAS (LOCATION AND VELOCITY OF MOLECULES) 

NADR1 = 0 
NADR2 = 0 
DO 200 I = 1 , NRD 
DO 200 J = 1 , NAD 
NM1 = IC(1,I,J) 

DO 205 N = 1 , NM1 
NADR1 = NADR1 + 1 
CALL RANDU(PP) 

P1(4,NADR1) = C(19,I, J)+DRX(PP-.5) 

CALL RANDU(P) 

PI ( 5, NADR1 ) = C(20,I, J)+DDALFAx(P-.5) 

DO 205 NV = 1,3 
203 CALL RANDU(P) 

V = -3.+6.XP 
B = EXP(-VXV) 

CALL RANDU(P) 

IF(B.LT.P) GO TO 203 
PI ( NV, NADR1 ) = PXSIN(B)XVTER1 
IFCNV. EQ.1)P1(NV,NADR1)=P1(NV,NADR1)+V0X 
I F( NV . EQ . 2) PI ( NV, NADR1 ) = P1 ( NV, NADRD+VOY 

205 CONTINUE 

REPEAT PROCEDURE FOR SPECIES 2 
NM2 =IC( 2, I , J ) 

DO 210 N = 1 , NM2 
NADR2 = NADR2+1 
CALL RANDU(P) 

P2( 4, NADR2 ) = C( 19 , I , J )+DRX( P- . 5) 

CALL RANDU(P) 

P2(5, NADR2) = C( 20 , I , J )+DDAL FAx( p- . 5 ) 

DO 210 NV = 1,3 
207 CALL RANDU(P) 

V = -3.+6.XP 
B = EXP(-VXV) 

CALL RANDU(P) 

IF(B.LT.P) GO TO 207 
P2 ( NV , NADR2 ) = PXSI N ( B ) XVTER2 
IF(NV.EQ.1)P2(NV,NADR2)=P2(NV,NADR2)+V0X 
IFCNV. EQ.2)P2(NV,NADR2)=P2(NV,NADR2)+V0Y 

210 CONTINUE 

WHEN NECESSARY, REPEAT PROCEDURE FOR SPECIES 3. 

200 CONTINUE 



SIM02890 
SIM02900 
SIM02910 
SIM02920 
SIM02930 
SIM02940 
SIM02950 
SIM02960 
SIM02970 
SIM02980 
SIM02990 
SIM03000 
SIM03010 
SIM03020 
SIM03030 
SIM03040 
SIM03050 
SIM03060 
SIM0307 0 
SIM03080 
SIM03090 
SIM03100 
SIM03110 
SIM03120 
SIM03130 
SIM03140 
SIM03150 
SIM0316 0 
SIM03170 
SIM03180 
SIM03190 
SIM03200 
SIM0321 0 
SIM03220 
SIM03230 
SIMO 3240 
SIM03250 
SIM03260 
SIM03270 
SIM03280 
SIM03290 
SIM03300 
SIM03310 
SIM03320 
SIM03330 
SIM03340 
SIM03350 
SIM03360 
SIM0337 0 
SIM03380 
SIM03390 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM03400 



C 

C 

C 



C 



C 



C 



DEFINE HERE ALL COLLISION PARAMETERS TO BE INCLUDED IN SIMULATION SIM03410 



RESET COLLISION TIMERS 
DO 19 1=1, NRD 
DO 19 J = 1 , NAD 
DO 19 L = 1 ,3 
DO 1 9M=1 , 3 
KT=3x(L-1)+M+9 
C(KT,I,J)=0. 

SET EXPECTED MAXIMUM RELATIVE VELOCITY 
KV=3X(L-1)+M 
EM1=SPEC( L , 1 ) 

EM2=SPEC(M, 1) 

EMR=EMlxEM2/( EM1+EM2) 

IF( KS . EQ . 1 )GO TO 17 
TEM=F0E1(4,I,KR,KS) 

GO TO 18 

17 TEM=FWP1(4,I,KR) 



BLOCK 9 SIM03420 
SIM03430 
SIM03440 
SIM03450 
SIM0346 0 
SIM03470 
SIM03480 
SIM03490 

IN COLLISIONS SIM03500 

SIM03510 
SIM03520 
SIM03530 
SIM03540 
SIM03550 
SIM0356 0 
SIM03570 
SIM03580 
SIM03590 
SIM036 00 
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18 RV=2./SQRTCPI*2.*B0LTZ*TEM/EMR) 

19 C(KV,I, J)=2.XRV 



SIM03610 
SIM03620 

C SIM03630 
C MAXIMUM RELATIVE VELOCITY WILL BE RESET IF FASTER ENCOUNTERS OCCUR SIM03640 
C SIM03650 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM03660 



TIME=0. 

LOOP OVER TIME INTERVALS BLOCK 

DO 6000 JDTM = 1,NIS 



SIM0367 0 
10 SIM03680 
SIM036 90 

C SIM03700 

C SIM03710 

TIME=TIME+DTM SIM03720 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM03730 



C 

C 

C 

C 



MOVE ALL MOLECULES BLOCK 

MOVE MOLECULES OF SPECIES 1 

DO 310 II = 1 , NMOL 1 



SIM03740 
12 SIM03750 
SIM037 6 0 
SIM0377 0 
SIM03780 



C SKIP INACTIVE MOLECULES 13 SIM53790 



IF CP1C4,I1) .EQ.-99. ) GO TO 310 
VX = Pl(l, II) 

V Y = PI C 2 , II ) 

RX = PICA, ID 
T = P1C5, II ) 

TOLD=T 



SIM03800 
SIM038 1 0 
SIM03820 
SIM03830 
SIM038A0 
SIM03850 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM03860 



1 A SIM0387 0 
SIM0388 0 
SIM03890 
SIM03900 
SIM0391 0 
SIM03920 
SIM03930 
SIM039A0 
SIM03950 
SIM0396 0 
SIM0397 0 
SIM0398 0 
SIM03990 
SIMOAOOO 
SIM0A010 
SIM0A020 
SIM0A030 
SIMOAOAO 
SIM0A050 
SIM0A06 0 
SIM0A07 0 
SIM0A080 
SIM0A090 
SIM0A100 
SIM0A1 10 
SIM0A120 
SIM0A130 
SIMOA1AO 
5IM0A150 
SIM0A160 
SIM0A17 0 
SIM0A180 
SIMOA1 90 
SIM0A200 
SIMOA21 0 
SIM0A220 
SIM0A230 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0A2A0 

C DEACTIVATE ALL MOLECULES THAT MOVED OUT BLOCK 17 SIM0A250 

301 IFCRNEW.LT. RMAX . AND . RNEW . GT . RMIN . AND . T . LT . TMAX . AND . T . GT . TMIN) GO TOSIM0A260 
X303 SIM0A27 0 

RI = C RNEW-RMIN )/DR+ .999999 SIM0A28 0 

IR=RI SIM0A290 

TI=CTMAX-T)/DDALFA+. 999999 SIM0A300 

IT =TI SIM0A31 0 

C SIM0A320 



C FIND NEW COORDINATES BLOCK 

XX=RXXCOSCT)+VXXDTM 

YY=RXXSINCT)+VYXDTM 

RNEW=SQRTCXXXX2+YYXX2) 

T=ATANC YY/XX) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C FOR LAST SECTOR FIND COLLISIONS WITH THE WALL AND SAMPLE THEM 

C BLOCK 15,16 

IFC ALFA . GT . DALFA)G0 TO 301 
IFCT.GT.O . )GO TO 301 
DTR=DTM*T/CT-TOLD) 

C DTR IS THE TIME REMAINING AFTER A MOLECULE STRIKES THE WALL 
I FC DTR . LT . 1 E-10 ) DTR = 1 E-10 
RW=RX+C RNEW-RX )XDTR/DTM 
IFCRW.LT. RMIN) RW=RMIN+DTRx .001 
IFCRW. GT. RMAX) RW=RMAX-DTR* . 001 
LOC=RW/DR+l . 

C 

C COUNT COLLISIONS WITH THE WALLCMUST BE WAIGHTED =XDFI) 

SSC 1 , LOC, JSAMP ) =SSC 1 , LOC, JSAMP)+1 

C 

C SET VELOCITY AND LOCATION AFTER A MOLECULE STRIKES THE WALL 
302 CALL RAN DU CP) 

B = VWM1XSQRT C-ALOGCP)) 

CALL RANDUCP) 

BB=2 . XPI XP 

P1C1,NADR1)=BXC0SCBB) 

VX=BXCOSCBB) 

P1C2,NADR1)=B*SINCBB) 

VY=B*SINCBB) 

CALL RANDUCP) 

P1C3,NADR1) =VWM1XSQRT C-ALOGCP)) 

XX=RXXC0SCT)+VXXDTR 

YY=RX*SINCT)+VYXDTR 

RNEW=SQRTCXXXX2+YYXX2) 

T=ATANC YY/XX) 
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IF(IT.LE.0)G0 TO 304 
IF( IT . GT . NAD )G0 TO 305 
I F( IR . GT . NRD )G0 TO 306 

C COUNT S DIRECTION, SAMPLE 

FSN1C1,IT,KR)=FSN1(1,IT,KR)+1. 

GO TO 309 

C COUNT W DIRECTION, SAMPLE 

304 IFCIR. LE. 0)IR=1 

I F( IR . GT . NRD) IR=NRD 
FWN1 C 1 , IR) =FWN1 ( 1 , IR)+1 . 

GO TO 309 

C COUNT E DIRECTION, SAMPLE 

305 IFCIR. LE.0)IR=1 
IFCIR. GT.NRD)IR=NRD 
FEN1C1,IR)=FEN1C1,IR)+1. 

GO TO 309 

C COUNT N DIRECTION, SAMPLE 

306 FNN1C1,IT,KR)=FNN1(1,IT,KR)+1 . 

SET NEW VALUES IN THE MOLECULE TABLE 



.18 

.18 



18 



.18 



309 RNEW=-99. 

303 PI C 4 , 1 1 ) =RNEW 
PI C 5 , 1 1 ) =T 

310 CONTINUE 



REPEAT PROCEDURE FOR SPECIES 2 MOLECULES. 

DO 320 12=1 , NM0L2 
SKIP INACTIVE MOLECULES , 

IFCP2C4, 12) . EQ .-99 . )GO TO 320 
VX = P2 C 1 , 12 ) 

VY = P2C 2, 12) 

RX = P2C 4 , 12 ) 

T = P2C5, 12) 

TOL D=T 



.BLOCKS 12 T018 



XX=RX*COSCT)+VX*DTM 
YY=RX*SINCT)+VYXDTM 
RNEW=SQRT CXXXX2+YYXX2) 
T=ATANC YY/XX) 
COLLISIONS WITH THE WALL 



DTR 



IFC ALFA . GT . DAL FA )G0 TO 311 
IFCT.GT.O. )GO TO 311 
DTR=DTM*T/CT-TOLD) 

IS THE TIME REMAINING AFTER A MOLECULE STRIKES THE WALL 
IFCDTR.LT. 1E-10)DTR=IE-10 
RW = RX+C RNEW-RX) XDTR/DTM 
IFCRW.LT. RMIN) RW=RW+DTRX .001 
IFC RW . GT . RMAX) RW = RW-DTR* .0 01 
LOC=RW/DR+l 

COUNT COLLISIONS WITH THE WALL CMUST BE WEIGHTED =*DFI) 

SSC 2, LOC, JSAMP )=SSC2,L0C, JSAMP) + 1 



FIND THE NEW COORDINATES OF THE MOLECULE 
312 CALL RANDUCP) 

B=VWM2*SQRTC-AL0G(P)) 

CALL RANDUCP) 

BB=2 . *PI*P 

P2C1,NADR2)=B*C0SCBB) 

VX=B*COSCBB) 

P2C2,NADR2)=B*SINCBB) 

VY=B*SINCBB) 

CALL RANDUCP) 

XX=RX*COSCT)+VX*DTR 

YY=RX*SINCT)+VY*DTR 

RNEW=SQRTCXX**2+YY*X2) 

T=ATANCYY/XX) 



DEACTIVATE MOLECULES THAT MOVED OUT. COUNT FOR OUTPUT FLUX EVALUATION 
311 IFCRNEW.LT. RMAX . AND . RNEW . GT . RMIN . AND . T . LT . TMAX .AND.T.GT.TMIN)GO 
*313 



SIM04330 
S IM04340 
SIM04350 
SIM0436 0 
SIM0437 0 
SIM04380 
SIM04390 
SIM04400 
SIM04410 
SIM04420 
SIM04430 
SIM04440 
SIM04450 
SIM0446 0 
SIM04470 
SIM04480 
SIM04490 
SIM04500 
SIM0451 0 
SIM04520 
SIM04530 
SIM04540 
SIM04550 
SIM04560 
SIM0457 0 
SIM04580 
SIM04590 
SIM04600 
SIM046 10 
SIM04620 
SIM04630 
SIM04640 
SIM04650 
SIM04660 
SIM0467 0 
SIM04680 
SIM04690 
SIM047 00 
SIM04710 
SIMO4720 
SIM04730 
SIM04740 
SIM047 50 
SIM04760 
SIM04770 
SIMO4780 
SIM047 9 0 
SIM04800 
SIM04810 
SIM04820 
SIM04830 
SIM04840 
SIM04850 
SIM04860 
SIM0487 0 
SIM04880 
SIM04890 
SIM04900 
SIM049 1 0 
SIM04920 
SIM049 30 
SIM04940 
SIM04950 
SIM0496 0 
SIM04970 
SIM04980 
SIM04990 
SIM05000 
SIM05010 
SIM05020 
TOSIM05030 
SIM05040 
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RI = CRNEW-RMIN)/DR+. 999999 SIM05050 

IR = RI SIM0506 0 

TI = (TMAX-T)/DDALFA+. 999999 SIM05070 

IT = TI SIM05080 

IF (IT.LE.O)GO TO 314 SIM05090 

IF (IT.GT.NAD)GO TO 315 SIM05100 

I F( IR . GT . NRD)GO TO 316 SIM05110 

C COUNT S DIRECTION, SAMPLE 18 SIM05120 

FSN2( 1 ,IT,KR)=FSN2C1,IT,KR)+1. SIM05130 

GO TO 319 SIM05140 

C COUNT W DIRECTION, SAMPLE 18 SIM05150 

314 IF( IR . LE. 0 )IR=1 S1M05160 

IF(IR.GT.NRD)IR=NRD SIM05170 

FWN2C 1 , IR) =FWN2( 1 , IR)+1 . SIM05180 

GO TO 319 SIM05190 

C COUNT E DIRECTION, SAMPLE 18 SIM05200 

315 IFCIR. LE. 0)IR=1 SIM05210 

I F( IR . GT . NRD) IR=NRD SIM05220 

FEN2( 1 , IR) =FEN2( 1 , IR)+1 . SIM05230 

GO TO 319 SIM05240 

C COUNT N DIRECTION, SAMPLE, 18 SIM05250 

316 FNN2C 1,IT,KR)=FNN2C1,IT,KR)+1. SIM05260 

C SIM0527 0 

319 RNEW = -99. SIM05280 

C 19 SIM05290 

313 P2C4,I2)= RNEW SIM05300 

P2(5, I2)= T SIM0531 0 

320 CONTINUE SIM05320 

C SIM05330 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM05340 



C NOW NEW MOLECULES WILL BE INTRODUCED 

C TOTAL JET FLUX INTO THE REGION WAS DETERMINED BY F1CI) MOL. /SEC. 

C OR FWP1 C I , KR) , FWP2( I , KR ) FOR EACH SPECIES 

C NUMBER OF MOLECULES TO BE ACTIVATED PER DTM IS F(I)XDTM PER CELL. 

C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

DO 390 1=1, NRD 
C INPUT 'W' MOLECULES 

IFCKS.EQ.DGO TO 322 
ANEW1=F0E1(1,I,KR,KS-1) 

ANEW2=F0E2C1,I,KR,KS-1) 

GO TO 323 

322 ANEW1=FWP1 ( 1 , I , KRlXDTM . 

ANEW2=FWP2( 1,1, KR )XDTM 

323 NEW1 =ANEW1 
NEW2=ANEW2 
REM1=ANEW1-NEW1 
CALL RANDUCP) 

IFCREM1 .GT.P)NEW1=NEW1+1 
P.EM2 = ANEW2-NEW2 
CALL RANDUCP) 

IFC REM2 . GT . P )NEW2 = NEW2+1 
C ACTIVATE NEW INPUT MOLECULES 

C TIME, LOCATION AND VELOCITY COMP. OF NEW MOLS. ARE RANDOM FUNCTIONS 

C SPECIES 1 

IFCNEWl.LT. l)GO TO 341 
C 

DO 340 11=1, NEW1 
CALL RANDUCP) 

ATIME=PXDTM 
CALL RANDUCP) 

RSTART=CFLOATCI-l)+P)XDR+RMIN 
VELX=FWP1C2, I , KR) 

VELY=FWP1 C 3, 1 , KR) 

TEM=FWP1 C 4 , 1 , KR) 

C THERMAL VELOCITY 

VTER1=S0RTC2.XB0LTZXTEM/SPECC1,1)) 

C TEMPERATURE IS TRANSFERRED THROUGH FWP1C4,I,KR) 
XX=RSTARTXCOSCALFA)+VELXXATIME 
YY=RSTARTXSINCALFA)+VELYXATIME 
RX=SQRTCXXxx2+YYXX2) 

T=ATANCYY/XX) 



SIM05350 
SIM05360 
SIM0537 0 
SIM05380 
SIM05390 
SIM05400 
SIM0541 0 
SIM05420 
SIM05430 
SIM05440 
SIM05450 
SIM05460 
SIM0547 0 
SIM05480 
SIM05490 
SIM05500 
SIM05510 
SIM05520 
SIM05530 
SIM05540 
SIM05550 
SIM0556 0 
SIM0557 0 
SIM05580 
SIM05590 
SIM056 00 
SIM056 1 0 
SIM05620 
SIM05630 
SIM05640 
SIM05650 
SIM05660 
SIM0567 0 
SIM05680 
SIM056 90 
SIM057 00 
SIM057 1 0 
SIM05720 
SIM057 30 
SIM05740 
SIM057 50 
SIM05760 
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IF( RX . LT . RMIN . OR . RX . GT . RMAX)G0 TO 340 
IF(T . LT . TMIN . OR . T . GT . TMAX)G0 TO 340 
C FOR MORE ACCURATE CALCULATION SET THESE MOLECULES IN OUTPUT FLOWS 
C 

C DEFINE THE VELOCITY COMPONENTS 
CALL RANDU(P) 

B=VTER1*SQRTC-AL0GCP)) 

CALL RANDU(P) 

BB=2 . XPIXP 
VELX=VELX+B*COSCBB) 

VELY=VELY+B*SINCBB) 

CALL RANDU(P) 

VELZ=VTER1XSQRT C-ALOGCP)) 

C DEFINE THE NEW MOLECULE TO BE ACTIVATED 
DO 325 IACT=1 , NM0L1 
IFCP1C4,IACT) .EQ.-99. )GO TO 326 

325 CONTINUE 

IF THERE IS NO ROOM FOR ADDITIONAL MOLECULE PRINT 'ALARM' 

IF( I ACT . GE . NMOL1 )GO TO 3004 

326 PI C 4 , IACT) =RX 
PI C 5 , IACT)=T 
P1C1,IACT)=VELX 
PI ( 2, IACT) =VEL Y 
PI ( 3, IACT)=VELZ 

340 CONTINUE 



REPEAT PROCEDURE FOR SPECIES 2 
341 I F( NEW2 . LT . 1 )GO TO 351 
DO 350 12=1, NEW2 
CALL RANDU(P) 

ATIME=PXDTM 
CALL RANDU(P) 

RSTART=RMIN+( FLOAT CI-1)+P)*DR 
VELX=FWP2C2, I,KR) 

' VELY=FWP2C3,I,KR) 

TEM = FWP2C 4 , 1 , KR ) 

THERMAL VELOCITY 

VTER2=SQRTC2.XB0LTZ*TEM/SPECC2,1)) 
XX=RSTARTXCOS(ALFA)+VELXXATIME 
YY = RSTARTXSINC ALFA )+VELYXATIME 
RX=SQRTCXXX*2+YY*X2) 

T=ATAN(YY/XX) 

IFCRX.LT. RMIN. OR. RX .GT . RMAX)GO TO 350 
IFCT.LT.TMIN.OR.T.GT. TMAX)GO TO 350 
FOR MORE ACCURATE RESULTS SET THESE MOLECULES IN OUTPUT FLOWS 

DEFINE VELOCITY COMPONENTS 
CALL RANDUCP) 

B=VTER2*SQRT C-ALOGCP)) 

CALL RANDUCP) 

BB=2 . xpixp 
VELX=VELX+BXC0SCBB) 

VELY=VELY+B*SINCBB) 

CALL RANDUCP) 

VELZ=VTER2*SQRTC-AL0GCP)) 

FIND A NEW MOLECULE TO BE ACTIVATED 
DO 345 IACT=1 , NM0L2 
IF CP2C4,IACT) .EQ.-99. )GO TO 346 

345 CONTINUE 

IF THERE IS NO PLACE THEN PRINT ALARM 
IFCIACT.EQ.NM0L2)G0 TO 3004 

346 P2C 4 » IACT ) =RX 
P2C 5, IACT)=T 
P2C 1 , IACT ) =VELX 
P2C2,IACT)=VELY 
P2C 3 , IACT) =VELZ 

350 CONTINUE 

REPEAT PROCEDURE FOR SPEC-3 



INPUT E MOLECULES 



SIM0577 

SIM0578 

SIM0579 

SIM0580 

SIM0581 

SIM0582 

SIM0583 

SIM0584 

SIM0585 

SIM0586 

SIM0587 

SIM0588 

SIM0589 

SIM0590 

SIM0591 

SIM0592 

SIM0593 

SIM0594 

SIM0595 

SIM0596 

SIM0597 

SIM0598 

SIM0599 

SIM0600 

SIM0601 

SIM0602 

SIM0603 

SIM0604 

SIM0605 

SIM0606 

SIM0607 

SIM0608 

SIM0609 

SIM0610 

SIM0611 

SIM0612 

SIM0613 

SIM0614 

SIM0615 

SIM0616 

SIM0617 

SIM0618 

SIM0619 

SIM0620 

SIM0621 

SIM0622 

SIM0623 

SIM0624 

SIM0625 

SIM0626 

SIM0627 

SIM0628 

SIM0629 

SIM0630 

SIM0631 

SIM0632 

SIM0633 

SIM0634 

SIM0635 

SIM0636 

SIM0637 

SIM0638 

SIM0639 

SIM0640 

SIM0641 

SIM0642 

SIM0643 

SIM0644 

SIM0645 

SIM0646 

SIM0647 

SIM0648 
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\ 



351 IF CTMIN.LE. 0. )G0 TO 390 



ANEW1 =F0W1 ( 1 , I , KR, KS+1 ) 



NEW1=ANEW1 
REM=ANEW1-NEW1 
CALL RANDUCP) 
IFCP.LT.REM)NEW1=NEW1+1 

ANEW2=F0W2C1,I,KR,KS+1) 
NEW2=ANEW2 
REM=ANEW2-NEW2 
CALL RANDUCP) 

IFCP. LT . REM) NEW2=NEW2+1 



IFCNEWl . LE. 0)GO TO 370 

DO 362 11=1, NEW1 
CALL RANDUCP) 

T=TMIN+DDALFA*P 
CALL RANDUCP) 

RNEW=RMIN+CFL0ATCI-1)+P)XDR 
CALL RANDUCP) 

TEM=F0N1CA,I,KR,KS+1) 

VTER1=SQRTC2.*B0LTZ*TEM/SPECC1,1)) 

B=VTER1*SQRTC-AL0GCP)) 

CALL RANDUCP) 

BB=2.*PI*P 

VELX=F0W1 C 2, 1 , KR, KS+1 )+B*C0SCBB) 
VELY=F0W1C3,I,KR,KS+1)+B*SINCBB) 

CALL RANDUCP) 

VELZ=VTER1*SQRTC-AL0GCP)) 

FIND A NEW MOLECULE TO BE ACTIVATED 
DO 36 A I ACT=1 , NMOL1 
IFCP1CA,IACT).EQ.-99.)G0 TO 366 
36 A CONTINUE 

IF THERE IS NO ROOM FOR ADDITIONAL MOLECULE THEN PRINT 'ALARM' 
IFC IACT . EQ . NMOL 1 )GO TO 300A 
366 PICA, IACT)=RNEW 
PI C 5 , IACT) =T 
PI C 1 , IACT)=VELX 
PI C 2, IACT) =VELY 
362 PI C 3 , 1 ACT) =VELZ 
370 CONTINUE 

REPEAT FOR SPECIES 2 

IFC NEW2 . LE . 0 )GO TO 390 
DO 382 12=1, NEW2 
CALL RANDUCP) 

T = TMI N+DDAL FAXP 
CALL RANDUCP) 

RNEW=RMIN+CFLOATCI-l)+P)*DR 
CALL RANDUCP) 

TEM=F0W2CA, I,KR,KS+1) 

VTER2=SQRTC2.*B0LTZ*TEM/SPECC2, 1) ) 

B=VTER2*SQRTC-AL0G<P)) 

CALL RANDUCP) 

BB=2 . XPIXP 

VELX=F0W2C2,I,KR,KS+1)+B*C0SCBB) 

VELY=F0W2C3, 1 , KR, KS+1 )+B*COSC BB) 

CALL RANDUCP) 

VELZ=VTER2*SQRTC-AL0GCP)) 

C FIND A NEW MOLECULE TO BE ACTIVATED 
DO 38A IACT=1 , NM0L2 
IFCP2CA, IACT) .EQ.-99)G0 TO 386 
38A CONTINUE 

C IF THERE IS NO ROOM FOR ADDITIONAL MOLECULE THEN PRINT ALARM 
IFC IACT . EQ . NMOL 2) GO TO 300A 
386 P2C A, IACT)=RNEW 



SIM06A90 
SIM06500 
SIM0651 0 
SIM06520 
SIM06530 
SIM065A0 
SIM06550 
SIM06560 
SIM0657 0 
SIM06580 
SIM06590 
SIM066 00 
SIM066 1 0 
SIM06620 
SIM06630 
SIM066AO 
SIM06650 
SIM06660 
SIM06670 
SIM06680 
SIM06690 
SIM067 00 
SIM06710 
SIM06720 
SIM067 30 
SIM067A0 
SIM067 50 
SIM0676 0 
SIM0677 0 
SIM06780 
SIM06790 
SIM06800 
SIM06810 
SIM06820 
SIM06830 
SIM068A0 
SIM06850 
SIM0686 0 
SIM0687 0 
SIM06880 
SIM06890 
SIM06900 
SIM06910 
SIM06920 
SIM06930 
SIM069A0 
SIM06950 
SIM06960 
SIM06970 
SIM06980 
SIM06990 
SIM07000 
SIM07010 
SIM07020 
SIM07030 
SIM070A0 
SIM07050 
SIM07060 
SIM07070 
SIM07080 
SIM07090 
SIM07100 
SIM07110 
SIM07120 
SIM07130 
SIM071A0 
SIM07150 
SIM07160 
SIM07170 
SIM07180 
SIM07190 
SIM07200 
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P2(5,IACT)=T 
P2( 1 , I ACT ) =VELX 
P2(2, IACT)=VELY 
382 P2C 3, IACT) =VELZ 



REPEAT FOR SPEC-3 
390 CONTINUE 



INPUT *S' MOLECULES 
DO A 60 J=1,NAD 
IF( KR . LT . 2)G0 TO 460 
ANEW1=FNN1(1, J,KR-1) 

IF(ANEW1.LE. .00001) GO TO 420 
NEW1=ANEW1 
REM=ANEW1-NEW1 
CALL RANDU(P) 
IF(P.LT.REM)NEW1=NEW1+1 
DO 402 11=1, NEW1 
CALL RANDU(P) 

T=TMIN+DDALFA*(P+FL0ATCJ-1)) 

CALL RANDU(P) 

RNEW=RMIN+PXDR 
CALL RANDU(P) 

TEM=FNN1 ( 4, J , KR-1 ) 

VTER1=SQRT(2 .XB0LTZXTEM/SPECC1,!)) 
B=VTER1XSQRT(-AL0G(P)) 

CALL RANDU(P) 

BB=2.*PI*P 

VELX=FNN1(2, J,KR-1)+BXC0S(BB) 
VELY=FNN1(3, J,KR-1)+BXSIN(BB) 

CALL RANDU(P) 
VELZ=VTER1*SQRT(-AL0G(P)) 

C FIND A NEW MOLECULE TO BE ACTIVATED 
DO 404 I ACT = 1 , NMOL 1 
IF( PI ( 4 , IACT ) .EQ.-99. )G0 TO 406 
404 CONTINUE 

I F( IACT . EQ . NMOL 1 )G0 TO 3004 
406 Pl(4, IACT)=RNEW 
Pl(5, IACT)=T 
Pl(l, IACT)=VELX 
P1C2, IACT)=VELY 
402 PI (3, IACT)=VELZ 
420 CONTINUE 
C REPEAT FOR SPECIES 2 

ANEW2=FNN2( 1 , J , KR-1) 

IF(ANEW2.LT. .OOOODGO TO 440 

NEW2=ANEW2 

REM=ANEW2-NEW2 

CALL RANDU(P) 

IF( P . LT . REM) NEW2=NEW2+1 
DO 422 12=1, NEW2 
CALL RANDU(P) 

T = TMIN+DDAL FAX ( p+ FLOAT ( J-l ) ) 

CALL RANDU(P) 

RNEW=RMIN+PXDR 
CALL RANDU(P) 

TEM=FNN2( 4 , J , KR-1 ) 

VTER2 = SQRT(2.XB0LTZXTEM/SPEC(2, D) 
B=VTER2XSQRT(-AL0G(P)) 

CALL RANDU(P) 

BB=2.XPI*P 

VELX=FNN2(2, J,KR-1)+BXC0S(BB) 
VELY=FNN2(3, J,KR-1)+BXSIN(BB) 

CALL RANDU(P) 

VELZ=VTER2xSQRT ( -ALOG( P) ) 

C FIND A NEW MOLECULE TO BE ACTIVATED 
DO 424 IACT=1 , NM0L2 
IF( P2( 4, IACT ) . EQ . -99 )G0 TO 426 
424 CONTINUE 

IF(IACT.EQ.NM0L2)G0 TO 3004 
426 P2(4,IACT)=RNEW 



SIM0721C 
SIM0722C 
SIM0723C 
SIM0724C 
SIM0725C 
SIM0726C 
SIM0727C 
SIM0728C 
SIM0729C 
SIM0730C 
SIM0731C 
SIM0732C 
SIM0733C 
SIM0734C 
SIM0735C 
SIM0736C 
SIM0737C 
SIM0738C 
SIM0739C 
SIM07 40( 
SIM0741C 
SIM0742C 
SIM0743( 
SIM0744C 
SIM0745C 
SIM0746C 
SIM0747C 
SIM0748C 
SIM0749C 
SIM0750C 
SIM0751C 
SIM0752( 
SIM0753C 
SIM07 54( 
SIM0755C 
SIM07 5 6 ( 
SIM0757( 
SIM07 5 8 ( 
SIM0759I 
SIM0760( 
SIM076K 
SIM0762( 
SIM0763C 
SIM0764C 
SIM0765C 
SIM07661 
SIM0767C 
SIM0768C 
SIM0769C 
SIM0770C 
SIM0771C 
SIM0772C 
SIM0773C 
SIM0774C 
SIM0775C 
SIM0776C 
SIM0777 C 
SIM0778C 
SIM0779C 
SIM0780C 
SIM0781C 
SIM0782C 
SIM0783C 
SIM0784C 
SIM0785C 
SIM0786 C 
SIM0787C 
SIM0788C 
SIM0789C 
SIM07900 
SIM0791C 
SIM07920 



88 



no o ooooo 



P2C5,IACT)=T 
P2(1,IACT)=VELX 
P2(2,IACT)=VELY 
422 P2( 3 , I ACT) =VELZ 
440 CONTINUE 



REPEAT FOR SPEC-3 



INPUT ' N ' MOLECULES 

SPEC-1 

ANEW1=FSN1(1, J,KR+1) 

IFCANEWl.LE. .00001) GO TO 450 

NEW1=ANEW1 

REM=ANEW1-NEW1 

CALL RANDU(P) 

IF ( P . LT . REM )NEW1=NEW1+1 
DO 442 11=1, NEW1 
CALL RANDU(P) 

T =TMIN +DDALFAX( P+ FLOAT ( J-l ) ) 

CALL RANDU(P) 

RNEW=RMAX-PXDR 
CALL RANDU(P) 

TEM=FSN1(4,J,KR+1) 
VTER1=SQRT(2.XB0LTZXTEM/SPEC(1,1)) 
B=VTER1*SQRT(-AL0G(P) ) 

CALL RANDU(P) 

BB=2 . xpixp 

VELX=FSN1(2,J,KR+1 )+BXCOS(BB) 
VELY=FSN1(3, J,KR+l)+BxSIN(BB) 

CALL RANDU(P) 

VELZ=VTER1*SQRT C-ALOGCP) ) 

FIND A NEW MOLECULE TO BE ACTIVATED 
DO 444 IACT=1 , NM0L1 
IF(P1(4, I ACT) . EQ.-99 . )GO TO 446 
444 CONTINUE 

I F( I ACT . EQ . NMOL1 )GO TO 3004 
446 Pl(4, IACT)=RNEW 
PI (5, IACT)=T 
PI ( 1 , 1 ACT) =VELX 
PI ( 2, 1 ACT) =VELY 
442 P1(3,IACT)=VELZ 
450 CONTINUE 



INPUT N MOLECULES SPEC-2 
ANEW2=FSN2( 1 , J , KR+1 ) 

I FC ANEW2 .LE. .00001) GO TO 460 
NEW2=ANEW2 
REM=ANEW2-NEW2 
CALL RANDU(P) 

IF(P . LT . REM)NEW2=NEW2+1 
DO 452 12=1, NEW2 
CALL RANDU(P) 

T =TMIN+DDAL FAX ( P+ FLOAT ( J-l) ) 

CALL RANDU(P) 

RNEW=RMAX-PXDR 
CALL RANDU(P) 

TEM=FSN2(4, J,KR+1) 
VTER2=S0RT(2.XB0LTZXTEM/SPEC(2,1)) 
B=VTER2XS0RT(-AL0G(P) ) 

CALL RANDU(P) 

BB=2XPIXP 

VELX=FSN2(2, J,KR+1)+BXC0S(BB) 
VELY=FSN2(3, J , KR+1) +BXSINC BB) 

CALL RANDU(P) 
VELZ=VTER2XS0RT(-AL0G(P)) 

C FIND A NEW MOLECULE TO BE ACTIVATED 
DO 454 IACT=1 , NM0L2 
IF(P2(4,IACT) .EQ.-99)G0 TO 456 
454 CONTINUE 

I F ( I ACT . EQ . NMO L 2 ) GO TO 3004 
456 P2(4,IACT)=RNEW 
P2( 5 , IACT ) =T 
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P2( 1 , IACT ) =VELX 
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P2(2, IACT)=VELY 
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452 


P2(3,IACT)=VELZ 
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CONTINUE 
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REARRANGE MOLECULES IN THEIR CELLS 
INITIALIZATION 
5000 CONTINUE 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
KIP=NM0L1+NM0L2+NM0L3 
DO 1001 KADD=1 , KIP 
1001 IP(KADD)=0 
KADD=0 
N 1 A = 0 
N2A = 0 
N3A = 0 

N1A,N2A,N3A ARE THE NUMBER OF ACTIVE MOLECULES (COUNTED NEXT) 

DO 1100 1=1, NRD 
DO 1100 J =1 , NAD 

DO 1005 K=1 , 4 
1005 IC(K,I,J)=0 

SET SPECIES 1 
N1C = 0 
N2C=0 
N3C = 0 
NTC = 0 

DO 1020 K1 =1 , NM0L1 
IF(P1(4,K1) .EQ.-99 . )G0 TO 1020 
RL0C=ABS(C(19,I,J)-P1(4,K1)) 

TLOC=ABS(C(20,I, J)-P1(5,K1)) 

IF(RL0C.GT.DRX.5)G0 TO 1020 
IF( TLOC . GT . DDALFA* . 5 )G0 TO 1020 

N1C=N1C+1 
NT C = NTC+1 
N1A=N1A+1 

IF(NTC.EQ.1)IC(4,I,J)=KADD 
IC(1,I,J)=IC(1,I, J)+l 
KADD=KADD+1 
IP ( KADD ) =K1 

1020 CONTINUE 

SET SPEC-2 MOLECULES 

DO 1040 K2=l , NM0L2 
IF(P2(4,K2) .EQ.-99. )GO TO 1040 
RL0C=ABS(C(19,I, J)-P2(4,K2)) 

TLOC=ABS(C(20,I,J)-P2(5,K2)) 

IFCRL0C.GT.DRX.5)G0 TO 1040 
IF(TLOC. GT. DDALFAX . 5)G0 TO 1040 

N2A=N2A+1 

NTC=NTC+1 

IF(NTC.EQ.1)IC(4,I,J)=KADD 
IC(2,I,J)=IC(2,I,J)+1 
KADD=KADD+1 
IP(KADD)=K2 
1040 CONTINUE 
1100 CONTINUE 

WRITE(6,1021)N1A,N2A 

1021 FORMAT ( ' NUMBER OF ACTIVE MOLECULES SPEC1= M5,' SPEC2=',I5/) 



SIM087 

SIM087 

SIM087 

SIM087 

SIM087 

SIM087 

-SIM087 

SIM088 

SIM088 

SIM088 

SIM088 

S.IM088 

SIM088 

SIM088 

SIM088 

SIM088 

SIM088 

SIM089 

SIM089 

SIM089 

SIM089 

SIM089 

SIM089 

SIM089 

SIM089 

SIM089 

SIM089 

SIM090 

SIM090 

SIM090 

SIM090 

SIM090 

SIM090 

SIM090 

SIM090 

SIM090 

SIM090 

SIM091 

SIM091 

SIM091 

SIM091 

SIM091 

SIM091 

SIM091 

SIM091 

SIM091 

SIM091 

SIM092 

SIM092 

SIM092 

SIM092 

SIM092 

SIM092 

SIM092 

5IM092 

SIM092 

SIM092 

SIM093 

SIM093 



XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM093 



CALCULATE COLLISIONS BLOCK SIM093 
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DO 999 1=1, NRD 
DO 999 J=1 , NAD 
999 NCOL C I , J ) =0 
C 

DO 900 1=1, NRD 
DO 900 J=1 , NAD 
C 

NM( 1)=IC(1,I,J) 

NMC2)=ICC2,I, J) 

NM( 3)=IC(3,I,J) 

NT0T=NM(1)+NM(2)+NM(3) 

DO 900 L = 1 , 3 
DO 900 M=1 , 3 

C NUMBER OF SPECIES IN THE PROGRAM IS 2 (3) 

KV=CL-1)*3+M 

KT=KV+9 

C KV, KT ARE THE ADDRESSES OF RELATIVE VELOCITY AND COLLISION TIMERS 
920 IF(C(KT,I, J) .GE.TIME) GO TO 900 
C C(KT,I,J)IS THE INTEGRATED TIME FOR L-M COLLISION 
KSEL=0 
KREJ=0 

IF(NM(L).GT.l.AND.NM(M).GT.l)GO TO 912 
C NO COLLISIONS ARE CALCULATED IF THERE ARE NO MOLECULES 

911 C(KT,I,J)=C(KT,I»J)+DTM 
GO TO 900 

C SELECT NOW THE MOLECULES FOR COLLISION 

912 IF CKSEL.GE.100)GO TO 911 
CALL RANDUCP) 

MOLl=PXNM(L)+. 999999 
IF(MOL1.EQ.O)MOL1=1 

C 

CALL RANDUCP) 

M0L2=PXNM(M)+. 999999 
IF(MOL2.EQ.O)MOL2=1 
C 

KSEL=KSEL+1 

CHECK IF THE SAME MOLECULE HAS BEEN SELECTED TWICE 
IF(L.EQ.M.AND.M0L1.EQ.M0L2)G0 TO 912 



FIND THE ACTUAL ADDRESSES OF THE SELECTED MOLECULES 

IFCL . EQ . 1 )K1 =0 
IF(L.EQ.2)K1=NM(1) 

IF(L.EQ.3)K1=NM(1)+NM(2) 

I FCM . EQ . 1 )K2=0 
IF(M.EQ.2)K2=NM(1) 

I FCM. EQ.3)K2=NM(1)+NM(2) 

KAD1=M0L1 + K1 + ICC 4, 1 , J ) 

KAD2=M0L2+K2 + ICC A , I , J ) 

KADI , KAD2 ARE THE LOCATION OF SELECTED MOLECULES IN IPC ) 
MAD1=IPCKAD1) 

MAD2 = I P C KAD2 ) 

MAD1 , MAD2 ARE THE ACTUAL ADDRESSES OF THE SELECTED MOLECULES CTHE 
INDICATION OF WHAT SPECIES THEY ARE HAS BEEN DEFINED HERE) 

DO 930 N=1 , 3 

IFCL . EQ . 1 ) VN1 =P1 C N, MAD1 ) 

IFCL .EQ.2)VN1=P2CN,MAD1) 

IFCL. EQ.3)VN1=P3CN,MAD1) 

IFCM. EQ.1)VN2=P1CN, MAD2 ) 

IFCM.EQ.2)VN2=P2CN,MAD2) 

IFCM.EQ.3)VN2=P3CN,MAD2) 

930 VRCCN)=VN1-VN2 
C VRCC3) CONTAIN THE THREE RELATIVE VELOCITY COMPONENTS 
VR = SQRT C VRCC 1 )XVRCC 1 )+VRCC 2 )XVRCC 2 )+VRCC 3) XVRCC3 ) ) 

C VR IS THE RELATIVE SPEED IN A SPECIFIC COLLISION 
IFCCCKV,I,J).LT.VR)CCKV,I,J)=VR 

C LAST STATEMENT RESETS THE MAXIMUM RELATIVE VELOCITY FOR FURTHER 
C CALCULATIONS 
C 

IFCKREJ.GT.100)GO TO 911 
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CALL RANDUCP) 

AVR=VR/CCKV,I, J) 

KREJ=KREJ+1 
IFCAVR.LT. P)G0 TO 912 

LAST STATEMENT REJECTS THE CALCULATED COLLISION 
NOW A SPECIES L-M COLLISION HAS BEEN SELECTED 

CALCULATE NOW THE PROBABILITY THAT SUCH A COLLISION WILL BE COUNTED 
FOR THE L AND M SPECIES RESPECTIVELY 
LP = 1 
LM=1 

CALL RANDU(P) 

ANM= FLOAT CNMCL))/FLOATCNMCM)) 

IFCANM.GT. 1 . )G0 TO 950 
IFCANM.LT. P)MP=0 
GO TO 955 



950 ANM=1./ANM 

IFCANM.LT. P)LP=0 



955 CXS=PI*CSPECCL,2)+SPECCM,2))x*2/4. 

ALP=LP 

ALM=LM 

V0LUME=REGC4,KR,KS)XV0LCI) 

CHECK 

DNL=FLOATCNMCL))/ VOLUME 
DNM=FLOATCNMCM) )/V0LUME 
USE EQ.10.3 

T0CCL,M)=ALP/CCXSXDNMXVRXNMCL))+ALM/CCXSXDNLXVRXNMCM)) 

SET THIS VALUE INTO CCKT,I,J) 

CCKT,I, J)=CCKT,I, J)+T0CCL,M) 



SAMPLE THIS COLLISION 

NCOLCI, J)=NC0LCI, J)+l 

FIND RELATIVE MASSES 
CALL RANDUCP) 

BB=1 .-2.XP 
AA=SQRTC1 .-BBxBB) 

VRCC 1 ) =BBxVR 
CALL RANDUCP) 
BB=2.XPIXP 

VRCC2)=AAXC0SCBB)XVR 

VRCC3)=AAXSINCBB)XVR 

FIND RELATIVE MASSES 

SM=SPECCL,1)+SPECCM,1) 

RML=SPECCL,1)/SM 

RMM=SPECCM,1)/SM 
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CALCULATE HERE THE ACTUAL ADDRESSES OF COLLIDING MOLECULES AND SET 
THEIR NEW VELOCITIES CNEW VELOCITY COMPONENTS ARE ADDED TO THE 
VELOCITY OF THE CENTER OF MASS VCCM) 

DO 960 N=1 , 3 
IFCL.EQ.1)V1=P1CN,MAD1) 

IFCL.EQ.2)V1=P2CN»MAD1) 

IFCL. EQ.3)V1=P3CN,MAD1) 

IFCM. EQ.1)V2=P1CN,MAD2) 

IFCM.EQ.2)V2=P2CN,MAD2) 

IFCM.EQ.3)V2=P3CN,MAD2) 

VCCM=RMLXV1+RMMXV2 
VCCM1 =VCCM+VRCC N ) xRMM 
VCCM2= VCCM -VRCC N ) XRML 
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C CHANGE IN VELOCITY IS INPUTED ONLY IF PROBABILITY OF COLL ISION . GT . 1 
IFCLP. NE.DGO TO 961 
IFCL.EQ.1)P1CN,MAD1)=VCCM1 

IFCL. EQ.2)P2CN,MAD1)=VCCM1 
I FC L . EQ . 3 ) P3 C N, MAD1) =VCCM1 

961 IF(LM. NE.DGO TO 960 

IFCM. EQ.1)P1CN,MAD2)=VCCM2 
IFCM.EQ.2)P2CN,MAD2)=VCCM2 
IFCM.EQ.3)P3CN> MAD2) =VCCM2 

960 CONTINUE 
GO TO 920 
900 CONTINUE 

WRITE! 6, 4545)T0CC1, 1 ) , T0CC 1 , 2) , TOC! 2, 1) ,T0CC2,2) DIME 
4545 FORMAT C * TIME* , 5E10 . 3) 
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xxxxxxxxxxxxxxxxxxEND OF COL LISIONS************x******xxx*xx**xx*xx*x**siMl 0980 
NON NEW TEMPERATURES MAY BE CALCULATED IN EACH CELL 
AVERAGE VELOCITY IN CELLS 



WRITE! 6 , 4547 ) 

9547 FORMAT! '0 I J 
1 VAVY1 VAVX2 



N1E 



N2E 

VAVY2 



TEMPR1 



TEMPR2 



VAVX1 



NCOL 1 ,/) 
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DO 1110 1=1, NRD 
DO 1110 J=1 , NAD 
N1E=IC( 1 , I , J ) 

N2E=IC(2, I , J ) 

KAD1=IC( 4 , I , J ) 

KAD2=ICC4,I, J)+N1E 
IFCN1E.LT. l)GO TO 1112 

VX1=0. 

VX2= 0 . 

VY1 =0 . 

VY2=0 . 

DO 1111 IM=1,N1E 
MAD1 =KAD1+IM 
IAD=IPCMAD1 ) 

VX1=VX1+P1(1,IAD) 

VY1 =VY1+P1 ( 2 , IAD ) 

IFCN2E. LE. 0)GO TO 1110 

DO 1115 IM=1 , N2E 

MAD2=MAD1+IM 

IAD=IPCMAD2) 

VX2=VX2+P2( 1 , IAD) 

VY2=VY2+P2( 2, IAD) 

AVERAGE VELOCITIES IN THE CELL WILL RESULT- 

VAVX1 =VX1/N1 E 

VAVY1 =VY1/N1 E 

VAVX2=VX2/N2E 

VAVY2=VY2/N2E 

IFCI. NE.DGO TO 1116 

FSN1C2, J,KR)=VAVX1+FSN1(2, J,KR) 

FSN1C3, J,KR)=VAVY1+FSN1C3, J,KR) 

FSN2C 2 , J , KR ) =VAVX2 + FSN2( 2 , J , KR) 

FSN2C 3, J , KR ) =VAVY2 + FSN2( 3 , J , KR) 

IFCI . NE.NRD)GO TO 1117 
FNN1C2, J,KR)=VAVX1+FNN1C2, J,KR) 

FNN1C3, J ,KR)=VAVY1+FNN1C3, J,KR) 

FNN2C 2, J , KR ) =VAVX2 + FNN2C 2 , J , KR) 

FNN2C 3, J , KR)=VAVY2+FNN2C 3 , J , KR) 

IFCJ. NE.DGO TO 1118 
FEN1 C 2 , 1 )=VAVX1 + FEN1 C 2 , 1 ) 
FEN1C3,I)=VAVY1+FEN1C3,I) 
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FEN2C2,I) =VAVXZ+FENZ( 2,1) 
FEN2C3,I)=VAVY2+FEN2C3,I) 

C 

1118 IFCJ.NE.NAD)GO TO 1119 
FWN1C2,I)=VAVX1+FWN1C2,I) 
FWN1C3,I)=VAVY1+FWN1C3,I) 
FWN2C2,I)=VAVX2+FWN2C2,I) 

FWN2C3, I)=VAVY2+FWN2C3,I) 

1119 CONTINUE 

THERMAL VELOCITIES AND TEMPERATURES 
ENRG1=0 . 

ENRG2=0 . 

IFCN1E.LT. 1)G0 T01131 
DO 1130 IM=1,N1E 
MAD1 =KAD1+IM 
IAD=IP(MAD1 ) 

CX1=P1C1, IAD1-VAVX1 
CY1=P1C2, IAD1-VAVY1 
CZ1=P1 C 3 , I AD) 

1130 ENRG1=ENRG1+CX1XCX1+CY1XCY1+CZ1*CZ1 

1131 IFCN2E.LT. 1)G0 TO 1150 
DO 11 AO IM=1,N2E 
MAD2=KAD2+IM 
IAD=IPCMAD2) 

CX2=P2C1,IAD)-VAVX2 
CY2=P2C2,IAD)-VAVY2 
CZ2=P2C 3, IAD) 



11 AO ENRG2=ENRG2+CX2XCX2+CY2XCY2+CZ2XCZ2 
C 

1150 TEMPR1=CENRG1/N1E)XSPECC1,1)/C3.XB0LTZ) 
TEMPR2 = C ENRG2/N2E) XSPECC 2,l)/C3.xB0LTZ) 

C 

IFCI.NE.DGO TO 1151 

FSN1CA, J,KR)=TEMPR1+FSN1CA, J,KR) 

FSN2CA, J,KR)=TEMPR2+FSN2CA, J,KR) 

~1151 IFC I . NE . NRD)GO TO 1152 

FNN1CA, J,KR)=TEMPR1+FNN1CA, J,KR) 

FNN2CA, J,KR)=TEMPR2+FNN2CA, J,KR) 

C 

1152 IFCJ.NE.DGO TO 1153 

FWN1 C A , I ) =TEMPR1+FWN1 C A , I ) 

FWN2C A , I ) =TEMPR2+FWN2C A , I ) 

r* 

“1153 IFC J . NE . NAD)GO TO 115A 

FEN1CA, I)=TEMPR1+FEN1CA,I) 

FEN2CA, I)=TEMPR2+FEN2C A, I) 

1 15A CONTINUE 
C 

WRITEC6,A5A8)I, J , N1E, N2E, TEMPR1 , TEMPR2, VAVX1 
1LCI, J) 

A5A8 FORMAT ( * • , AI5, 6 F12 . 3, 17 ) 

C 

1110 CONTINUE 

n 

^6000 CONTINUE 

C CALCULATE AVERAGED PARAMETERS, WEIGHTED BY DFI 
DO 6010 1=1, NRD 
DO 6010 KPAR=1 , A 

FOW1CKPAR, I,KR,KS)=FWN1CKPAR,I)/CNIS) 

F0W2C KPAR, I , KR, KS ) = FWN2C KPAR, I )/C NIS ) 

FOE1 C KPAR, I , KR, KS ) =FEN1 C KPAR, I)/CNIS) 

6010 F0E2CKPAR,I,KR,KS)=FEN2CKPAR,I)/CNIS) 

C 

DO 6020 J=1 , NAD 
DO 6020 KPAR=1 , A 

FSN1 C KPAR, J , KR ) = FSN1 CKPAR,J,KR)/NIS 
FSN2CKPAR, J , KR ) =FSN2 C KPAR, J,KR)/NIS 
FNN1CKPAR, J,KR)=FNN1CKPAR, J,KR)/NIS 



6020 FNN2(KPAR, J , KR)=FNN2(KPAR, J,KR)/NIS SIM12250 

C SIM12260 

C TO PREPARE FLOWS FOR THE NEXT REGION OR SECTOR, DIVIDE FONO BY L0CALSIM12270 
C DFI AND MULTIPLY BY NEXT DFI WHEN STARTING NEW REGION SIM12280 

C SIM12290 

C CALCULATE HERE THE MEAN FREE PATH AND STORE INTO REG( , , ) SIM12300 

C SIM1231 0 

C STOP PROGRAM IF FLOW BECOMES COLLISIONLESS OR NUMBER DENSITY IS EQUAL SIM12320 
C TO THE AMBIENT NUMBER DENSITY. STORE F0E1,F0E2 IN A SEPARATE FILE TO SIM12330 

C BE USED IN A DIFFERENT PROGRAM. SIM12340 

C SIM12350 

C PRINT AVERAGED RESULTS SIM12360 

WRITE( 6 , 6 029 ) SIM12370 

6029 FORMAT ( ' ’,////) SIM12380 

WRITE(6,6030)NIS,KR SIM12390 

6030 FORMAT ( ' AVERAGED OUTPUT FLOWS AFTER’, 15,' TIME INCREMENTS IN REGISIM12400 

ION’, 15) SIM1241 0 

WRITE(6 ,6031) SIM12420 

6031 FORMAT ( ’ 0 I FEN1C1,I) FEN2(1,I) FWN1C1,I) FWN2( 1 , I ) ’ SIM12430 



1) 

DO 6033 1=1, NRD 

WRITEC 6 ,6032)1, FEN1 C 1,I),FEN2(1,I), FWN1 ( 1,I),FWN2(1,I) 

6032 FORMAT ( ' ’,I5,4F13.3) 

6033 CONTINUE 
WRITE (6,6034) 

6034 FORMAT ( * 0 J FNN1(1,J,KR) FNN2(1,J,KR) FSN1(1,J,KR) 

12(1, J,KR) ') 

DO 6036 J=1 , NAD 



SIR12440 
SIM12450 
SIM12460 
SIM1247 0 
SIM12480 
SIM12490 
FSNSIM12500 
SIM12510 
SIM12520 



WRITE(6,6037)J,FNN1(1,J,KR),FNN2(1,J,KR),FSN1(1,J,KR),FSN2(1,J,KR)SIM12530 



6037 FORMAT ( ’ ’,I5,4F15.5) 
6036 CONTINUE 



SIM12540 
SIM12550 

C SIM1256 0 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX SI Ml 257 0 

SIM12580 
SIM12590 
SIM12600 
SIM12610 
SIM12620 
SIM12630 



C START A NEW REGION 
C IF(KR.EQ.IO) GO TO 7000 

C KR=KR+1 

C GO TO 2000 

C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM12640 



C SIM12650 

7000 CONTINUE SIM12660 

C SIM12670 

C SIM12680 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX SIM12690 
C SIM12700 

C FIND IF FLOW BECAME COLLISIONLESS IN THE WHOLE SECTOR SIM12710 

C IF POSITIVE, STORE F0E1,F0E2 IN A SEPARATE FILE AND STOP PROGRAM SIM12720 

C SIM12730 

C PREPARE DATA FOR THE NEXT SECTOR SIM12740 

C KR = 1 SIM127 50 

C STOP PROGRAM IF KS WAS BOUNDED BY THE WALL SIM12760 

C KS=KS+1 SIM12770 

C GO TO 1000 SIM12780 

C SIM12790 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM12800 



C SIM12810 
C IF THERE IS BACK FLOW (FWN1,FWN2) NONZERO CALCULATE NEXT ITERATION SIM12820 
C ITER=ITER+1 SIM12830 
C KS=KR=1 SIM12840 
C GO TO 3000 SIM12850 
C SIM1286 0 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM1287 0 



WRITE( 6 , 6 099 ) 

6099 FORMAT ( ’ 1 DATA FOR TEN MOLECULES SPEC. 2’) 

DO 3003 1=1,10 

WRITE( 6 , 3002 ) P2( 1 , 1 ) , P2( 2 , 1 ) , P2( 3, 1 ) , P2( 4 , 1 ) , P2( 5, 1 ) 

3002 FORMAT ( ' ’ , 4F1 0 . 3 , E18 . 1 0 ) 

3003 CONTINUE 
GO TO 3009 

3004 WRITE(6 ,3005) 



SIM12880 
SIM1289 0 
SIM12900 
SIM129 1 0 
SIM12920 
SIM12930 
SIM12940 
SIM12950 
SIM12960 
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oo 



3005 FORMAT ( * NO PLACE FOR ADDITIONAL 
3009 CONTINUE 
STOP 
END 



SUBROUTINE RANDU(P) 
COMMON IX 
IY = IXX65539 
IF (IY) 5,6,6 

5 IY = IY+21 47483647+1 

6 P = IY 

P = PX.4656613E-9 

IX = IY 

RETURN 

END 



MOLECULES*) SIM12970 

SIM12980 
SIM12990 
SIM13000 
SIM13010 
SIM13020 
SIM13030 
SIM13040 
SIM13050 
SIMI3060 
SIM13070 
SIM13080 
SIM13090 
SIM131 00 
SIM13110 
SIM13120 
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B.5 Program SIMUL - User's Guide 

Preparation: Run AXSYM program. From its output data evaluate: 

ALFA - The averaged angle for the continuum breakdown 

(P - 0.05) 

TETA - Flow direction along the breakdown limit. 

Flow parameters along this line - Pressure, temperature, velocity, 



free path. 

Input Data: 

Radius of nozzle ring R1 

Radial size of a cell DR 

Maximum radius in simulation RP 

Angle of breakdown limit ALFA 

Flow direction along (ALFA) TETA 

Averaged flow velocity along (ALFA) Vo 

Molecular weight of each species Spec(I,l) 

Molecular diameter of each species Spec (1,2) 

Mean free path along (ALFA) FPM 

Averaged Temperature (ALFA) TEMP 

Time increment DTM 

Number of time increments NIS 

Options 

a. Geometry 



The program is designed to run for axisymme trie ring flow. 



mean 



For a two dimensional planar flow the molecular motion, collisions and 
flow calculation remain unchanged. The flow cross section remain 



unchanged. Instead of the angle DFI the two dimensional flow requires 
the definition of the width of each region (or cell). To keep the 



number of molecules within reasonable computational limits this size 
has to decrease in the same manner as DFI. 

The part of the program which has to be changed for this purpose is 
lines 230 to 250 . 
b. Wall flux calculations: 

The program may run for the whole molecular region resulting the flux 
towards the wall (this part needs additional debugging). However in 
order to make it more efficient we may stop the program at a sector 
where the flow becomes collisionless. The remaining part of the flow 
may be regarded either as collisionless or if we define a much larger 
size of cells we may calculate the molecular collisions on this 
basis . 

This part need additional analysis and programming. 

Execution Commands 

Without additional changes the program runs under WATFIV compiler using 
the following command: 

WATFI AXSYM * (XTYPE) 

Further developments will be required to run the program for each sector 
separately and the output intermediate results on the mass storage. 
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B.6 



The Influence of the Ambient Gas 



The temperature, pressure and density of the ambient gas are shown in 
Table 1. Because its temperature is much higher than in the jet gas, the 
thermal velocity of ambient gas molecules is much higher than jet molecules, 
the following contribution may be expected due to the ambient gas: 

a. Collisions between the "hot" ambient molecules with the "cold" jet 
molecules may cause an increase in the dissipation in the outer 
layer of the jet and increase in the flux towards the walls. 

b. For higher ambient pressures, those collisions become rare because 
of the low number density therefore, the influence of the 
collisions may decrease. 

The only way to evaluate the influence of these two controversial factors 
is by an additional simulation program to be designed for this region. The 
following are the main factors to be included in this program: 

Boundary Conditions: 

a. The jet side: F0E1 and F0E2 obtained from the last simulated sector 

(SI1IUL) supply the number flux, flow velocity componetns and 
temperatures. These parameters are given for all points along the 
radius R(I) at constant distances DR. In the low density domain the 
resolution DR is much too large compared with the expected aiean free 
path. A different mesh has to be designed for this purpose. It is 
possible that one cell may be sufficient. 

b. Far Field Condition: The boundary conditions where the gas may be 

regarded "undisturbed" by the jet are as follows 
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- Jet gas molecules are allowed to go out the simulation region 
(these molecules will be regarded as "lost” molecules). 

- Ambient gas is allowed to enter the simulated region according 
with their thermal velocity and number density. 

c. Solid Wall Bondaries . The solid wall may be assumed to have a 

constant temperature T w . Incident molecules of either species are 
reflected back from the wall. Different models of collisions with 
the wall may be employed. 

- Elastic collisions f specular reglection* (this calculation has been 
included in SIMUL program) 

- Collisions with ideal heat transfer (diffuse reflections) 

- Other models depending on the materials and surface parameters. 

The collision with the wall was included as an internal routine in the 

program SIMUL. If the general molecular simulation contains the program 
proposed here the collision with the wall will be omitted from SIMUL and 
become the core of the additional collisionless program. Figure 25 shows the 
low density (collisionless) region and its boundaries. 
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Figure 25. The Low Density Region in the Jet 
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SUMMARY OF REPORT 



Algorithms for the continuum regime and for the region where molecular 
collisions are significant have been developed. 

Program AXSYM contains the calculation of planar jet flow and 
axisymmetric ring jet flow. This program supplies data for the limits where 
the continuum approach become invalid and molecular approach should be 
employed . 

Program SIMUL contains the molecular simulation for the axisymmetric ring 
flow. This program may run for the whole molecular region to result in the 
calculation of flux towards the solid wall. For a more efficient simulation 
it is proposed to design an additional program for the collisionless region 
where ambient gas may be included. 

For the two dimensional flow, program SIMUL may be used after changing 
the definition of the cell geometry. 

To run the whole program it may be required to make separate runs for 
each sector and store results on the mass storage. 
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